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A COMPUTER PROGRAM FOR THE GEOMETRICALLY 


NONLINEAR STATIC AND DYNAMIC ANALYSIS OF ARBITRARILY 
LOADED SHELLS OF REVOLUTION, THEORY AND USERS MANUAL 
by Robert E. Ball 

SUMMARY 


A digital computer program known as SATANS - Static and Transient 
Analysis, Nonlinear, Shells, for the geometrically nonlinear static and 
dynamic response of arbitrarily loaded shells of revolution is presented. 
Instructions for the preparation of the input data cards and other 
information necessary for the operation of the program are described in 
detail and two sample problems are included. The governing partial 
differential equations are based upon Sanders 1 nonlinear thin shell 
theory for the conditions of small strains and moderately small rotations. 
The governing equations are reduced to uncoupled sets of four linear, 
second order, partial differential equations in the meridional and time 
coordinates by expanding the dependent cariables in a Fourier sine or 
cosine series in the circumferential coordinate and treating the nonlinear 
modal coupling terms as pseudo loads. The derivatives with respect to 
the meridional coordinate are approximated by central finite differences, 
and the displacement accelerations are approximated by the implicit 
Houbolt backward difference scheme with a constant time interval. At 
every load step or time step each set of difference equations is repeatedly 
solved, using an elimination method, until all solutions have converged. 

All geometric and material properties of the shell are axisymmetric, but 
may vary along the shell meridian. The applied load may consist of any 
combination of pressure loads, temperature distributions and initial 
conditions that are symmetric about a datum meridional plane. The shell 
material is isotropic, but the elastic modulus may vary through the 
thickness. The boundaries of the shell may be closed, free, fixed, or 
elastically restrained. The program is coded in the FORTRAN IV language 
and is dimensioned to allow a maximum of 10 arbitrary Fourier harmonics* 
and a maximum product of the total number of meridional stations and the 
total number of Fourier harmonics of 200. The program requires 155,000 
bytes of core storage. 



INTRODUCTION 


The design of many shell structures is influenced by the geometrically 
nonlinear response of the shell when subjected to static and/or dynamic 
loads. As a consequence, a number of investigations have been devoted 
to the study of the buckling phenomenon exhibited by shells. Most of 
the early works examine the behavior of the shallow spherical cap, the 
truncated cone, and the cylinder under axisymmetric loads. As a con- 
sequence of the lack of information on the axisymmetric response of 
shells with other meridional geometries and on the response of shells 
subjected to asymmetric loads, a computer program for the geometrically 
nonlinear static and dynamic response of arbitrarily loaded shells of 
revolution has been developed. The dynamic analysis capability is a 
recent extension of the program developed by the author for the non^ 
linear static analysis of arbitrarily loaded shells of revolution^- 1 ]. 

The program can be used to analyze any shell of revolution for which the 
following conditions hold: 

1) The geometric and material properties of the shell are axi- 
symmetric, but may vary along the shell meridian. 

2) The applied pressure and temperature distributions and initial 
conditions are symmetric about a datum meridional plane. 

3) The shell material is isotropic, but the modulus of elasticity 
may vary through the thickness. Poisson’s ratio is constant. 

4) The boundaries of the shell may be closed, free, fixed, or 
elastically restrained. 

The governing partial differential equations are based upon Sanders’ 
nonlinear thin shell theory for the condition of small strains and 
moderately small rotations!- 2 ]. The inplane and normal inertial forces 
are accounted for, but the rotary inertial terms are neglected. The 
set of governing nonlinear partial differential equations is reduced 
to an infinite number of sets of four second-order differential 
equations in the meridional and time coordinates by expanding all 
dependent variables in a sine or cosine series in terms of the cir- 
cumferential coordinate. The sets are uncoupled by utilizing appro- 
priate trigonometric identities and by treating the nonlinear coupling 
terms as pseudo loads. The meridional derivatives are replaced by 
the conventional central finite difference approximations, and the 
displacement accelerations are approximated by the implicit Houbolt 
backward differencing scheme>-3J # This leads to sets of algebraic 
equations in terms of the dependent variables and the Fourier index. 

At each load or time step, an estimate of the solution is obtained 
by extrapolation from the solutions at the previous load or time steps. 

The sets of algebraic equations are repeatedly solved using Potters *L4] 
form of Gaussian elimination, and the pseudo loads are recomputed, 
until the solution converges. 
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An automatic variable load incrementing routine is included in the 
program for the static analysis. When the number of iterations are 
small the load is incremented in equal steps. As the nonlinear terms 
become large, and the number of iterations exceeds a prescribed maximum, 
the incremental load is reduced by a factor of five. Any number of 
increment reductions can be made. The load is continually increased 
until either the prescribed maximum number of load steps or increment 
reductions have been taken. Post-buckling behavior cannot be determined 
in the static analysis because of the method of solution employed. 

This report contains a description of the theory, the method of 
solution, instructions for the preparation of the input data cards, 
and other information necessary for the operation of the program. Two 
sample problems are included to illustrate the data preparation and 
output format. For additional information concerning the accuracy and 
applicability of the program refer to references [5] and [ 6 ], 
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SYMBOImD 


a « 

B 

b 

D « 

d 
E 
E 

o 

s’ Q’ s0 


s 


h 

h 

c 

K 


W k se 


reference length 

inplane stiffness, Jfedg/(l-v ) 

nondimensional inplane stiffness, B/ ( E 0 h Q ) 

"bending stiffness, J£ 2 Edg/(l-\^) 

Q 

nondimensional “bending stiffness, D/ (E Q h^) 
elastic modulus 
reference elastic modulus 

= nondimensional Fourier coefficients for the reference 
surface strains , equations ( 32 ) 

nondimensional Fourier coefficient for the transverse 
force , equations ( 32 ) 

nondimensional Fourier coefficient for the effective 

transverse force, § /(cr h ) 

s o o 

thickness 

reference thickness 

last meridian station on the shell 

= nondimensional Fourier coefficients for the bending 
strains , equations ( 32 ) 


M , M q , M q = bending and twisting moments per unit length 
s b sb 

m = mass density of the shell material 


W m se 


= nondimensional Fourier coefficients for bending and 
twisting moments, equations ( 32 ) 


“T 


nondimensional Fourier coefficient for the thermal 
bending moment, equation (32) 


N ,N~,N n = membrane forces per unit length 
S U So 

N p = effective shear force, equation (l4) 
s b 

n = Fourier index 
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p,p s ,p e 


= nondimens ional Fourier coefficients for the components 
of the pressure load, equations (32) 


Q s ,Q 0 

A. 

^S 

q s ?q e 


R ,R 
s’ 0 


r 


s 

T 

T 

( 

t 





U,v 


u ? v 


W 


- transverse forces per unit length 

= effective transverse force , equation (15b) 

,q = meridional 5 circumferential, and normal components of 
applied pressure load 

= principal radii of curvature 

= normal distance from the axis of the shell 

= meridional shell coordinate 

= time 

= reference time 

= nondimens ional time, T/T 

' o 

,t = nondimens ional Fourier coefficients for membrane 
S forces, equations (32) 

= nondimensional Fourier coefficient for the effective 

shear force, ft n /(cr h ) 
s0' o o 

= nondimensional Fourier coefficient for the thermal membrane 
force, equations (32) 

- displacements tangent to the meridian and to the parallel 
circle respectively. 

= nondimensional Fourier coefficients for the displacements 
tangent to the meridian and to the parallel circle 
respectively, equations (32) 

= displacement normal to the reference surface 


■w = nondimensional Fourier coefficient for the displacement 

normal to the reference surface, equations (32) 

a = coefficient of thermal expansion 


(3 = nondimensional coefficients for the nonlinear terms 

s in the strain-displacement relations, equations 
(20a) and (20c) 
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Y 


p'/p 


A 

6t 


nondimensional distance between stations, meridian 
length/ a/ (K^-l) 

nondimensional time interval 


V 


ss 

’es 


e0 5 e s0 = reference surface strains, equations (5) 

= thermal membrane force, equation (12c) 

= coordinate normal to the reference surface 

,TL,T 1 = nondimensional coefficients for the nonlinear terms in 
,T]g0 ? T| s 0 the equilibrium equations, equations (20c) 

- circumferential angle 


W H s0 


V = 

§ 


<u s ,<B e = 

§ S ’V § 





= bending strains, equations (6) 
thermal bending moment, equation (l2d) 
nondimensional mass, a jmd^/(h Q E o T o ) 

Poisson^ ratio 

nondimensional meridional coordinate, s/a 
nondimensional radius, r/a 
reference stress level 
nondimensional curvatures, a/P , a/p 

S u 

= reference surface rotations, equations (7) 

= nondimensional Fourier coefficients for the rotations, 
equations (32) 

local temperature change 


_ 3 _ 

30 

3 w 3 

9 s ° r 3 ? 

Fourier series coefficient 
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MATRICES 


A,B,B,C 

E,F,G,H,J 

e,f 

gjg 

A 

P ? P S P 

x 

z 


A,n 


A * £"2 


M- 


4x4 matrices, equations (25b) and (28b) 

4x4 matrices, reference [l] 

1x4 column matrices, reference [l] 

1x4 column matrices, equations (25b) and (28b) 

1x4 boundary condition matrix 
4x4 matrices, equations (30) 

1x4 column matrix, equations (29a) and (31a) 

1x4 column matrix containing the unknown variables 

u,v,w, and m 
s 

4x4 nondimens ional boundary condition matrices 
4x4 boundary condition matrices 
mass matrix 
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THEORY 


Shell Geometry 


Consider the general shell of revolution shown in figure 1. 
Located within this shell is a reference surface. All material points 
of the shell can be located using the orthogonal coordinate system s, 
9, £, where s is the meridional distance along the reference surface 
measured from one boundary, 0 is the circumferential angle measured 
from a datum meridian plane, and Q is the normal distance from the 
reference surface . The positive direction of each coordinate is indi- 
cated in figure 1. For convenience, let the reference surface be 
positioned so that 


J CEdC = 0 


( 1 ) 


where E is the elastic modulus and the integration is carried out over 
h, the thickness of the shell. Thus, when E is independent of £ the 
reference surface coincides with the middle surface of the shell. Further, 
let the location of the reference surface be described by the dependent 
variable r, the normal distance from the axis of the shell. Accordingly, 
the principal radii of curvature of the reference surface are 


R 0 = r/[l - (r') 2 ]= 

R = - [1 - (r') 2 ]2/ r " 

S 


( 2 ) 


where a prime denotes differentiation with respect to s. Further, note 
the Codazzi identity 


<1/ ■ r ' < E s 1 - E 0 1) / r 
and the relation 


( 3 ) 


r 


F! 


r/R R. 

1 S 0 


( 4 ) 


Strain-displacement Relations 


For a shell of revolution, the strain-displacement relations 
derived by Sanders take the form 


e = U’ + W/R + (§ 2 + § 2 )/2 
s s s 

eg = V*/r + r' u/r + W/Rq + (f 2 + § 2 )/2 


9 


s s0 = ^ V ' + U */ r " r ' V / r + § s V/ 2 


y ( 5 ) 


and 
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I 

I 


reference 

surface 



Figure 1. Shell Geometry and Coordinates 
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I 


H _ 5 ’ 

S "S 

* e = * q/ t + r ' § s / r 

"se * »a + ‘ r ' V r + iR t ' O i] / 2 


>• (S) 


where e , e~, and e are the reference surface strains, h - and 
s ? 0 ? s0 3 s 3 0 3 

K Q are the bending strains, U and V are displacements in the directions 
s 0 

tangent to the meridian and to the parallel circle respectively, W is 
the displacement normal to the reference surface, and § , § , and $ 
are rotations defined by y 


g = - w* + U/R 
s 7 s 

l e = - W/r + V/R 0 
§ = (V* + r'V/r - U‘/r)/2 


(7) 


In these equations, and henceforth, a superscript dot denotes differentia- 
tion with respect to 0. The positive direction of each displacement and 
rotation variable is indicated in figure 2. 


Equations of Motion 

Converting Sanders 1 equilibrium equations to the equations of 
motion for a shell of revolution leads to 

(rH s )' + r 0 - r’N 0 + rC^/R, + (R* 1 - R 0 X ) m; 0 /2 = r( 

~ + r(§ s N s + § 6 N s0^ R s + Ci(W s + V ] ‘/ 2 

N 6 + (rN sQ } ' + r ' N s0 + r V R 0 + rC(R 0 1 ~ = r(JmdOd 2 V/ST 2 

- rq 0 + r(S 0 R 0 + N s0 )/R q - r [•(*„ + ^)]’/2 

(rQ g )' + Q 0 - rR s /R g - rNg/Rg = r ( JWOc^/BT 2 

" r( l + N s + r§ 0 N s0 ) ’ + ( -s W s0 + *0 V* 


m 


and 

(rM s )' + M* Q - r M Q - rQ g = 0 (9) 

VHe 1 ' +r ' M se- r «e = 0 (10) 
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when the effects of rotary inertia are neglected. In equations (8) - 

(10), m is the mass density of the shell material, T is time, q , q a , 

s 0 

and. q are the meridional, circumferential, and normal components of the 
applied pressure load, and are the transverse forces per unit 

length, N , and N are the membrane forces per unit length, and 

S 0 S 0 

M , and M Q are the bending and twisting moments per unit length, 
s 0 s 0 

Refer to figure 3 for the positive directions of the pressure components, 
forces, and moments. 



e T = J* a t E d C / (1 - v) (12c) 

n T = J C a t E d £ / (1 - v) (I2d) 

In equations (12c) and (l2d), t is the local temperature change and a 
is the coefficient of thermal expansion. 


12 



Boundary Conditions 


In Sanders' nonlinear theory, the conditions to prescribe on the 
edge of a shell of revolution are 


A 


N 

s 

or U 

W S0 0r v 


or ¥ 

M or $ 
s s 


Y (13) 


where ft ^ and 15 are the effective shear and transverse forces per 
S0 s 

unit length defined by 


\e - + E e 1 - E = 1 )M se /2 + (N s + V * /2 


\ V 1 f a H se 


(HO 

( 15 a) 


Using the equilibrium equation (9) to eliminate from equation (15a) 
leads to 


= [ ( rM s } + 2M S0 - r M el r • $ s N s - § e N s9 


(15b) 


Elastic restraints at the edge of a shell can be provided for by linearly 
relating the forces or moment to the appropriate displacements or 
rotation. Consequently, the boundary conditions may be given in the 
matrix form 



N 

s 


U 


CD 


V 

Q 

§ 

^s 

+ A 

w 


$ 

s 


M 

s 

— — 


where Q and A are 4x4 matrices and Z is a column matrix. The values of 
the elements of these matrices are determined by the conditions pre- 
scribed at the shell boundary. 
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METHOD OF SOLUTION 


Fourier Expansions 

The crux of the method used here to solve the nonlinear field 
equations is the elimination of the independent variable 0 by expand- 
ing all dependent variables into sine or cosine series In the circum- 
ferential direction. Only loading and initial conditions that are 
symmetric about a datum meridian plane will be considered. Thus, the 
variable can be expressed in the form* 

00 

*, - r 1 V n> =° s ne (17> 

° n=0 


where a is a reference stress level, E is a reference elastic modulus, 
o o ’ 


and the nondimensional series coefficient cp 


(*) « 


is a function of the 


independent variables s and T. Similar series expansions can be made 
for the remaining dependent variables. 


Modal Uncoupling 


In order to eliminate the independent variable 0 from the problem, 

and convert the partial differential equations to sets of uncoupled 

partial differential equations, the nonlinear terms are treated as 

known quantities or pseudo loads. Since every nonlinear term is the 

product of two Fourier series, each product can be reduced to a single 

trigonometric series wherein the coefficient is itself a series. For 

example, using equation (17) 5 ^ can be expressed as 

s 


co oo 


» „ (n) 


§ s = Z Z *r' V- cos m0 cos n0 

° m=0 n=0 


Since 


cos m0 cos n 0 = ^ [cos (m-n) 0 + cos (m+n) 0] 


(18) 


( 19 ) 


equation (l8) can be given in the form 


•^Theoretically, the complete Fourier series including both the sine and 
cosine expansions should be used because of the possibility of "odd" 
displacements occurring under "even" loads, i.e. a bifurcation phenomenon. 
This aspect is not considered here. 
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s E 

o 


r 1 p s n) cos ne 


(20a) 


n=0 


where 


P‘ n> =&- l [^ 0 (i "' ) ^» s |i -” l 3 


(20b) 


° i=0 


with 

0 for n - 0 1 for i ^ n 

T1 = , H = 

1 for n > 0 2 for i = n 

Similar series expressions can he derived for the other nonlinear terms 
in equations (5)* (8), (l4), and (15b). They are 


t 2 =^ 

s e e 

o 


* E 

o 


§ $ = — 
s 9 E 

o 


i 

n=0 

00 

l 

n=0 


Pg cos n0 


P cos n0 


L p s e sinn0 


n=l 


IN = a h 
S s o o 


Vse = CT o h o 


§ N = a h 
s o o 


I N>s = a h 
0 o o 


) T] cos n0 
L ; S S 


n=0 


1 Ties cos ne 


n=0 

00 


I \ sl 


in n0 


n=l 


I "e sl 


sin n0 


r (2oc) 


n=l 


15 


Ve = CT o h o 1 V sin n0 


n=l 


§ N = ct h 
s s 0 o o 


I \e slnne 


n-1 


where h is a reference thickness . 
o 

As a result of the trigonometric series expansions, there is one 
set of governing equations for each value of n considered; when only 
the linear terms are considered the sets are uncoupled. The presence 

(n) 

of the nonlinear terms couples the sets through terms like f3' ' as 

s 

given by equation (20b). However, by treating the nonlinear terms as 
known quantities and grouping them with the load terms, the sets of 
equations become uncoupled. 


Final Equations 

f 7 1 

Budiansky and Radkowski have shown that for the linear shell 
problem each set of Sanders 7 uncoupled field equations can be reduced 
to four second order differential equations provided Mg is replaced 

by the equality obtained from the constituitive relations (lid) and 
(He) 

M q = vM s + D (1-v 2 )h 0 - (1-v)^ (21) 


to prevent derivatives of W higher than two from appearing. The 
same procedure is used here. The four unknown dependent variables 

are the nondimens ional series coefficients u^ n \ v^, w^ n \ and m 

? 5 5 s 


corresponding to U, V, W, and M respectively. Three of the final 

s 

four equations are derived from the equations of motion (8) by applying 

the rotational equilibrium equations (9) and (10), the constituitive 

relations (11 ) and (2l), and the strain-displacement relations (5), 

(6), and (7). The fourth equation is derived from the meridional 

bending moment -curvature relationship given by (lid) with k and r 

s y 

expressed in terms of the displacements. 


A convenient representation of these four equations is the non- 
dimensional matrix form 


,(n) »" + F (n) „(n)' 


+ G 


(n) z (n) , e (n) + - j2 a (n) /a( .2 (22) 


where 
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» - 


u 


w 


(n) 

r (n) 

(n) 

(n) 


t is the nondimens ional time T/T , T is a reference time, and p, is the 
mass matrix given hy 


\i = [i 


10 0 0 
0 10 0 
0 0 10 
0 0 0 0 


The nondimens ional scalar mass is defined hy 


^ = 


2 

a 


h E T 
o o o 


/ tnd C 


where a is a reference length. Henceforth, the superscript n will he 
dropped for convenience. 

The E, F, G, and e in equation (22) are matrices defined in 
reference [l]. The elements of E, F, and G are identical with those 
given in reference [7] for the linear shell analysis, hut the e matrix 
as defined in reference [l] contains both the load and thermal terms 
and the nonlinear terms . 

The boundary conditions on z are obtained hy applying the consti- 
tutive relations (ll) and (21), and the strain-displacement relations 
(5) 3 (6), and (7) to equation (l6). This leads to the matrix equation 

QHz T + (A + QJ) z = - Q f (23) 

where Q and A are the nondimens ional forms of Q and A. Matrices H and J 
are identical with those given in reference [7] for the linear shell 
problem, and matrix f, as defined in reference [l], contains the thermal 
and nonlinear terms. In this formulation, Q, A, and & are not functions 
of n, and hence, the same set of boundary conditions applies for each 
value of n considered. An example of the modifications required to 
allow different values of £ for each mode is given in reference [5]. 
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Spatial Finite Difference Formulation 

Let the shell meridian he divided into K - 1 equal increments, and 
denote the end of each increment or station by the index i. Thus, 
i - 1 corresponds to the initial edge of the shell and i « K corres- 
ponds to the final edge. One fictitious station is introduced off each 
end of the shell at i - 0 and i - K + 1. 


Let the first and second derivatives of z at station i be 
approximated by 

(24a) 

(24b) 


! i = ( Vi - 


z i = K-4/. - 2z -i + 


i+1 


J i-r 


where A is the nondimens ional distance between stations. Substituting 
equations (24a) and (24b) into equation (22) leads to 


A. z. _ + B. z. + C. z. _. = g. + 2A [L. (3 z/dt ). 

l l+l li li-l & i ^i v J i 


(25a) 


where 

A. = 2E./A + F- 

l i / i 

B. = - 4E./A + 2 A G. 

1 l' 1 

C. = 2E . /A - F. 

l l' l 


% ■ 2 & e i 


> (25b) 


and i = 1, 2 . . . K to insure equilibrium over the total length of 
the shell . 

At the boundaries equation (23) must be satisfied. Thus, 
substituting equation (24a) into equation ( 23 ) leads to 


k n i H i z 2 + (a L J i + V Z 1 - k 1 H 1 Z 0 = - l f i ( 26a ) 

at the initial edge, and 

2A Q K H K Z K+1 + + V Z K " 2A = \ ^ 2 ^ 


at the final edge . 
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Timewise Differencing Scheme 

The inertial terms that appear in equations (25) can be approxi- 
mated by HoubolVs backward differencing scheme. Accordingly, 


^\, s ■ (2 *i, 3 - 5z i,d-i + s.j-2 - 


(27) 


where j denotes the time step and 6t is the nondimensional time inter- 
val. Thus, substituting equation (27) into equation (25) yields 


A . z . . + B . z . . + C . z . . . = g . . 

1 l+l , 3 l 1,3 l 1-1,3 1,3 


( 28 a) 


where 


IjAii, 


B. = B. 
1 1 


(6t)‘ 


7 — g + 

S i,J 


2Au. 


(6t) c 


y (28b) 


(-5 z. . n + 4 z. . 0 - z. ) 
i,J-l ijJ-2 i,j-3 


and i = 1, 2, . . . K. 


Solution by Elimination 

Eqautions (26) and (28) constitute a set of simultaneous algebraic 

equations in the unknowns z. . provided g. . , z. , z. . 

1 j J J ^ 5 J ™ 1 j J ™ ^ ^ ana 

z. . 0 are known. There is one such set for each value of n considered. 

i* J-3 

The equations can be arranged in the form shown in figure 4. Since 
these equations are tridiagonal in the matrix sense. Potters 1 form of 
Gaussian elimination can be used to solve for the z. .. In this 
method, recursion relationships of the form 


. = P. g. . - P. x. . 

1,0 i ijo i 1 - 1 , a 

( 29 a) 

. . = - P. Z. . + X. . 

i,j i i+i j a ijj 

( 29 b) 


are developed. A forward pass from the initial edge to the final edge 
computes the x. ., and a back substitution determines the z. . . The 


19 


1 

J 

Ojjf ) + h 1 

o h ( n) 
°1 H 1 

2A 






Z (N) 
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• 

r (N) 

°K 
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b k 

4 B) 


Z (N) 
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2A 

Vk N) + 

vP 

2A 


M 

^K+l 


Wf’j 


Figure 4. Matrix Equation for n = N 
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matrices P^, P^, and P^ are independent of the load and solution. 
Hence , they are computed only once. They are 


- [°L “i 1 B 1 + J ll + "l] 1 C n i S (I + A lJ 


P. 

1 


[B. - C. P. . ] 
i i l-l 


-1 


P. 

l 


P. A. 

l l 


S (30) 


P. = P. C. 

l ii 


The initial value of x is 

X 1 ■ [°1 S 'i 1 B 1 * J i} + 'll 1 K + ^ bs c 

and the value of z at station K + 1 is 

JT _2_ 

Z K+1 = [°FC 2S - P K-1 V - (n K J K + V P J 

+ °K S ( - P K-l X K + Vl 5 ' (Q K J K 




(31a) 


+ V X K “ 



(31b) 


Poles 

The equations (26a) and (26b) are applicable when the shell has 
edges. If the shell has a pole, r=0, and special "boundary” condi- 
tions are required to assure finite stresses and strains at the pole. 
These conditions are derived in reference [l]. 


SOLUTION PROCEDURE 

As a consequence of the selection of the Houbolt timewise 
differencing scheme, both static and dynamic analyses can be carried 
out using essentially the same set of equations and solution procedure. 
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Static Analysis 


For a static analysis, p,=0, and the applied load is increased 
monotonically. Thus, the index j denotes the load step. 

The procedure used to determine z for the monotonically increasing 
load is illustrated in figure 5 and described below: 

1) The matrices P^, P^, and are computed. 

2) A solution is obtained for a specified increment (DEL0AD) 
of each Fourier coefficient of the design load. All pseudo loads 
are taken as zero. 

3) The new solution is used to calculate the nonlinear terms, 

and a new value of the load vector g is obtained for each n. Additional 
values of n may be introduced by the nonlinear terms. 

4) A solution is obtained for the new value of g for each n and 
is compared with the previous solution. 

5) If the difference between two consecutive solutions, at any 
station and for any n, is greater than a specified percentage (EPS) 
of the maximum solution in that mode then step #3 is repeated. How- 
ever, if the number of iterations has exceeded a specified maximum 
(ITRMAX), the total load (A10AD) is reduced by one load increment, 
the increment (DEL0AD) is reduced by a factor of 5? and this new 
increment is added to the load. If a specified number of load 
reductions (LCHMAX) have been made, the program ends. 

6) If the two consecutive solutions have converged, another load 
increment is added, provided the number of load steps is less than a 
specified maximum (LSMAX). An estimate of the solution for this new 
load is made by linear extrapolation using the two preceeding converged 
solutions, and step #3 is repeated. 

Since the method of solution is based on a nonlinear pseudo 
load approach, the shell reacts equally, in a linear fashion, to any 
change in either the applied load or the pseudo load. Thus, failure 
of the solution to converge in any mode can be attributed to two 
types of nonlinear behavior. Both types are illustrated in figure 5 . 

The existance of a maximum or an inflection point on the softening 
load-deflection curve A represents a type of behavior for which a 
solution can be obtained only below the points of zero or nearly zero 
slope. On the other hand, the existence of a stiffening nonlinearity, 
as illustrated by curve B of figure 5> can also cause a convergence 
failure when ever the slope becomes too steep. Thus, in general, it 
is necessary to examine the load -displacement behavior of the shell 
in order to determine the cause of the convergence failure . 
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Figure 5* Typical Load-Displacement Curves 
from a Static Analysis 
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Dynamic Analysis 


The dynamic analysis proceeds in essentially the same manner as 
the static analysis. The only differences are due to the fact that; 
(l) the applied load is not raonotonically increased, hut instead is a 
function of the time step j; and (2) initial conditions on z and 
dz/9t are required to start the procedure. A brief description of the 
procedure used to obtain the response of the shell for a specified 
period of time and time increment (DEL0AD) is given below: 

, — ^ 

1) The matrices P^, P^, and P^ are computed. 

2) The solutions at j = 0, -1 and -2 are computed for each n 
from the specified initial conditions using the expressions 


z. q = initial condition on z supplied by user 
(3z/3t). n = initial condition on 3z/3t supplied by user 
z i,-l ■ z i,0 - 6t (3z / 3t) i,0 
z i,-2 * z i,0 - < Sz / 3t >i,0 


for i = 0,1,... K+l. An estimate of the solution at j=l is obtained 
for each n from 


z i,l - z i,0 + 6t < 3z / St >i,0 


for i = 0, 1, 2, . . . K+l. 


3) This new solution is used to calculate the nonlinear terms, 
and a new value of g is obtained for each n using the estimated non- 
linear terms and applied loads at j and the solution at j-1, j-2, and 

j-3. 


4) A solution is obtained for the new value of g for each n and 
is compared with the previous solution at 

5) If the difference between two consecutive solutions, at any 
station and for any n, is greater than a specified percentage (EPS ) of 
the maximum solution in that mode then step #3 is repeated. However, 
if the number of iterations has exceeded a specified maximum (ITRAMAX) 
the program ends. 

6) If the two consecutive solutions are sufficiently close, an 
estimate of the solution at j+1 is obtained by quadratic extrapolation 
from the solution at j , J-l, and j-2 , The preceeding solutions are up- 
dated, and step #3 is repeated for the new time step j-j+1, provided 
the number of time steps is less than a specified maximum (LSMAX) . 
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Two comments are in order here. First, the approximations used 
to obtain the solutions at j = -1 and -2 are not the ones suggested by 
Houbolt. Eoubolt's approximations require a change in the B matrix 
at the first time step. This_is time consuming since it necessitates 
the recomputation of the P^, P^, and P^ matrices, and does not appear 

to be worth the extra effort. Second, the time interval is usually 
so small no iteration is required since the difference between the 
estimated solution and computed solution is generally negligible. 
However, when the shell becomes dynamically unstable, the solution 
may not converge, even with iteration. Thus, the maximum number of 
iterations allowed should be small. 


COMPUTER PROGRAM 
Brief Description 

The program described in this report - SATANS - Static and Transient 
Analysis, Nonlinear, Shells.,- is a modified version of the program 
described in reference (1)* The revisions were made by personnel at the 
NASA Langley Research Center and by the original author. The main 
difference between the two versions is the addition of the capability 
for dynamic analysis. Another difference is in the manner in which core 
storage is allocated for the solution vector z. The solution vector is 
now handled as a two dimensional array instead of a three dimensional 
array, allowing the user the freedom of prescribing almost any combination 
of meridional and circumferential unknowns within the dimensions of the 
array. In the modified program up to 200 unknowns may be specified so 
that the product of the total number of meridional stations and the total 
number of Fourier harmonics must be less than 201. However, the maximum 
number of Fourier harmonics that can be considered is still 10. Any 
combination of harmonics may be used. For example, n = 5, 0, 22, and 91 
is allowed; there is no restriction on the order nor on the number. 

A change was also made in the test for convergence. The original 
program required two consecutive solutions to differ by less than a 
specified percentage of the latest solution. This test was made at 
every station, for every mode, except when the solution was less than 
-6 

10 . Experience with this routine showed it to be too restrictive. 

Consequently, it was replaced with the requirement that for each 
harmonic the difference between two consecutive solutions at each 
station must be less than a specified percentage of the maximum solution 
in that harmonic, considering all the stations, except when the solution 
-5 

is less than 10 . This new test for convergence appears to provide 

converged, accurate solutions in fewer iterations than the original 
scheme . 

The output subroutine was also modified in order to present the 
data in more compact form; the COMMON and DIMENSION statements were 
changed to allow the compilation of the program in any order; and 
several bugs were detected and eliminated. The operational parameters 
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of the program and the "boundary conditions are still read in on cards, 
but the geometry and mass of the shell, the inplane and bending 
stiffnesses, the pressure and thermal loads, and the initial 
conditions are introduced through user- prepared subroutines. The 
input and output data may be in either dimensional form or non- 
dimensional form, and no special tapes, discs, or routines are 
required for execution. However, a tape is required if the dynamic 
analysis is to be restarted. All of these changes have enlarged 
the program to the extent that it now requires a core space of 
approximately 150,000 bytes on an IBM 360/67 digital computer and 
can no longer be executed on a 32,000 word couputer. The compilation 
time using the FORTRAN IV Compiler, Level H, is slightly less 
than 2 minutes on the NPS IBM 360 / 67 . The static version of this 
program has been available from COSMIC as M 7 O-IOO 98 , LAR- 10736. 

The computer program has been used to solve a number of static 
and dynamic problems for both axisymmetric and asymmetric loads [5*6]. 
Two of these problems are presented here to illustrate the input and 
output features of the program. 


Nond imens i onal i z a t i on 

The input and output data may be in either dimensional or non- 
dimensional form. The dimensional parameters are E q , a reference 

elastic modulus, a reference stress, a, a reference length, and 

h , a reference thickness. The variables are made nondimens ional as 
o 3 

follows : 


p 

= r/a 

t (n) = H (n) /(o h ) 

§ 

= s/a 

f (n) - 4 (n) /(a h ) 

Y 

■ al A/ a 

m (n) = M (n) (a/a h 3 ) 

UD 

S 

- a /B s 

u (n) = u( n )( Eo / aCTo ) 

"e 

- a / R e 

cp (n) = § ^ n ^(E /a ) 
^s s v o' 0' 

b 

- 

e (n) = « (n) (E >„ ) 
s s ^ o' o' 

d 

- »/(Vo 3) 

k ( n ) = H ( n )(aE /a ) 

s s x o' 0 ' 

M- 

- I”«(a 2 A 0 E 0 T o 2 > 

4 n) - % (n) (*w 

t (n) 

t T 

' *T <n) A» 0 i> 0 ) 

4 n) ’ ,, T (n)(a A„ h 0 3 ) 


> (32) 
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etc . 


Similar expressions hold for 




Example of a Static Analysis 

The first problem is the static analysis of a clamped, shallow 
spherical cap of constant thickness and uniformly loaded over one- 

half of the shell from 0 - -90° to 0 = 90° , This problem was first 
considered by Famili and Archer [8]. The geometry of the spherical 
cap can be specified by the single nondimens ional parameter where 

X = 2[3(l-v 2 )]^(H/h)^ 


H is the rise of the shell and h is its thickness. The classical 
buckling pressure of a complete sphere is denoted by q Q , where 

i 

q o = 2Eh 2 /R s 2 /[3(l-v 2 )] 2 

For this analysis, 

X =6 

v = .3 

meridian length = 105 in. 

R = R = 1000 in. 
s 0 

Z! o 

E = 27.3 x 10 lb/inf 
h =1 in. 

£ 

B = 30.0 x 10 lb/in. 

£ 

D = 2.5 x 10 lb - in. 

q = -30 lb/in. 2 -tt/2 s; 0 s n/2 

q o = -33.1 lb/in. 2 

The reference parameters for nondimens ionalizat ion are taken as 

6 P 

E ^ 30 x 10 lb/in. 
o ' 

a Q = 1000 lb/in. 2 
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a = 1000 in. 


h =' 1 in. 
o 

Seven stations over the length of the meridian and four modes are 
used for the solution. For the purpose of illustration, only the 
first three Fourier harmonics of the applied load are used. Thus, 

= -15.0 lb/in. 2 

q^ = -19.1 lb/in. 2 

q^ = 6.37 lb/in. 2 

The boundary conditions are 

U = V = ¥ = <£ = 0 

s 

This problem took, 33 minutes of execution time on the NPS IBM 360/67 
Computer using the FORTRAN IV, Level H, Compiler with 0PT=2. 


Example of a Dynamic Analysis 

The second example is the dynamic analysis of a clamped truncated 
cone subjected to an impulsive loading which is uniform along the 
meridian and varies in a cosine distribution over one -half of the 
circumference. This problem is Sample Problem No. 3 in the series of 
sample problems suggested by the Lockheed Missiles and Space Co.. The 
initial conditions are 


¥ = 0 . 

d¥/dT = -4440.8 cos 0 (in. /sec] 
dW/dT = 0. 

The physical parameters are 
v - .286 

meridian length = 15.004 in. 

r . = 7.9^99 in. 

ram 

r = 10.2300 in. 

max 

R = CD 

s 

Rq = r/cos 0) 


0^0^ 2tt 
-tt/2 £ 0 £ tt/2 
tt/ 2 s 0 ^ 3 tt/2 
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w = sin” kC 10. 2300 - 7.9^99)/l5.O04] 

E = 3.52 x 10^ lb/ in. 2 
h - .543 in. 

c 

B = 2.0816 x 10 lb/ in. 

D = 5-114 x lO^ lb - in. 

The time step is 

AT = 2 x 10 sec 

The reference parameters for nondimens ionaliz at ion are taken as 

E = 3.52 x 10^ lb/ in. 2 
o 

CT q = 1000 lb/ in. 2 

a = 15.004 in. 

h = .543 in. 
o 

T q = IO .965 x 10”^ sec 

Thirty- one stations over the length of the meridian and four modes 
are used for the solution. The first four Fourier harmonics of the 
initial conditions are 

dW^°V dT = -4440.8/tt in. /sec 

dW^VdT = -4440.8/2 in./sec 

dW^'/dT = -2 x 4440.8/(3™) in./sec 

dW^4 dT = 2 x 444 o.8/(15tt) in./sec 

The boundary conditions are 

U = V= W= § = 0 

s 

This problem took approximately 8 minutes of execution time for 
750 time steps on the WPS IBM 360/67 computer using the FORTRAN 
IV, Level H, Compiler with 0PT=2. 
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Card 

Columns 

Format 

Item 

Input Lata 

Static 

Example 

Cards 

Dynamic 

Example 

Interpretation 

1 

2-72 

i8a4 

TITLE 



Problem description. 

2 

1-5 

15 


10 

108 

The problem number. 

2 

6-10 

15 

S^RD 

0 

1 

Set to: 

1. = 0 for static analysis; 

2. ^ 1 for dynamic analysis. 

2 

11-15 

15 

IM0DE 

1 

0 

Set to: 

1. > 0 if modal data are desired 

for each harmonic; 

2. £ 0 if modal data are not 

desired. 

2 

16-20 

15 

NDIMEN 

0 

0 

Set to: 

1. = 0 if dimensional form of out- 
put data is desired; 

2. ^ 1 if nondimensional form of 
output data is desired. 

2 

21-25 

15 

NTHMAX 

1 

2 

The summed solution will be 
printed at NTHMAX meridians, 
0 £ NTHMAX s6, 

2 

26-30 

15 

IFREQ 

1 

13 

The solution will be printed at 
meridional stations 1, IFREQ + 1, 
2 * IFREQ + 1 , . . . . , and the final 
station. 

2 

31-35 

15 

IFRINT 

1 

2 

Every IFRINTth converged solution 
will be printed. 

2 

36-40 

15 

IBCINL 

-1 

0 

Set to: 

1 < 0 if the shell has a pole at 
the first station; 

2. ^ 0 if the shell does not have 
a pole at the first station. 

2 

41-45 

15 

IBCFNL 

0 

0 

Set to: 

1. < 0 if the shell has a pole at 

the final station; 

2. s 0 if the shell does not have 

a pole at the final station. 

2 

46-50 

15 

KMAX 

7 

31 

The number of meridional stations . 


The product of KMAX and MAXM (the 
number of Fourier terms in the 
solution) must he less than 201. 
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2 51-55 15 

2 56-60 15 

2 61-65 15 


2 66-70 15 


2 71-75 15 


2 76-80 15 


3 1-12 E12.3 

3 13-24 E12 . 3 


MNMAX 3 4 


maxm 4 4 


1SMAX 99 


750 


LCHMAX 2 


0 


Number of Fourier terms used to 
describe the initial conditions, 
the pressure loads, and the 
thermal loads , MNMAX ^ MAXM . 

Number of Fourier terms in the 
solution, MAXM £ 10 and (KMAX)* 
(MAXM) £ 201. 

Static analysis 

Maximum number of load intensities 
to be considered. For a nonlinear 
analysis this number should be 
large. For a linear analysis set 
ISMAX = 1. 

Dynamic analysis 

Maximum number of time increments 

to be considered, LSMAX=T /AT. 

max' 

Static analysis 

Maximum number of load increment 
reductions. He commend value, 2-4. 

Dynamic analysi s 
LCHMAX = 0 


ITRMAX 5° 20 Maximum number of iterations at 

any load intensity or time step. 
Recommended value, 10-30. For a 
linear analysis set ITRMAX = 1. 

IC 0 Static analysis 

IC = 0 


NU 

SIG0 


1 


.3 .286 
1000. 1000. 


Dynamic analysis 
Set to: 

1. ^ 0 if shell at rest at t = 0, 

or if restarting solution at 
t>0 , i.e. ITAFE = 2 or 3; 

2. > 0 for non-zero initial conditions 

at t=0. 

Poisson's ratio, v. 

Reference stress level, a . When 
the data is to be input in dimen- 
sional form set SIG0 = 1.. 
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3 


25-36 E12.3 ELAST 


3 37-48 E12.3 TKN 


3 49-60 E12 . 3 CHAR 

3 61-72 E12.3 TEE0 


4 1-12 E12.3 DELOAD 


•3E8 .35?E7 Reference modulus of elasticity, 

E . When the data is to be 
o 

input in dimensional form set 
ELAST a 1.. 

1. .543 Reference thickness, h Q . When the 

data is to be input in dimensional 
form set TKN =1.. 

1000. 15.004 Characteristic shell dimension, a. 

When the data is to be input in 
dimensional form set CHAR =1.. 


0 . 


Static analysis 

= O _ 


Dynamic analysis 

IO. 965 E -5 Reference time T . 

o 

.2 Static analysis 

The load increment. DELOAD remains 
unchanged until the solution fails 
to converge in ITRMftX iterations. 
Then it is automatically reduced 
by a factor of 5 . A maximum of 
LCHMAK reductions will occur 
provided the number of load 
intensities considered is less 
than LSMAX. 


1.82396E-2 Dynamic analysis 

The nondimens ional time increment 
fit. 


4 

5 


13-24 

E12.3 

EPS 

.01 

.01 

The convergence criterion. 
Recommended value, .01. 

1-5 

15 

ITAPE 

0 


Static analysis 
ITAPE = 0 


Dynamic analysis 
The parameter for obtaining the 
data to restart the solution at 
t > 0: 

1. no read or write on tape, ITAPE 

= 0 ; 

2. write Z, Zj6, Z2, and Z3 after 
final time step, I TAPE = 1; 

3. read Z,Zj6,Z2, and Z3 before 
initial time step, ITAPE = 2; 

4. read Z,Z)6,Z2, and Z3 before 
initial time step, and write 
Z,Z^,Z2, and Z3 after final 
time step, ITAPE = 3. 
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6 


1-72 


6E12.3 


0 . 


0 . 


3.14159 


A list of NTHMAX circumfer 
ential coordinates 0 , in 
radians, where print-out 
of the solution is desired 
This card is omitted if 
NTHMAX = 0. 


For a static analysis, the execution for each case terminates and the 
program transfers to the first read statement when either the number 
of load intensities considered equals ISMAX or the number of iterations 
equals ITRMAX and the number of previous DEL0AD reductions equals 
LCHMAX. 


For a dynamic analysis, the execution terminates when either 
ISMAX time increments have been taken or when the solution does not 
converge after ITRMAX iterations . 

Restart option. When Z, Z 0, Z2 and Z3 have been put on tape unit 
#8 after the final time step (ITAFE = l), the response computation can 
be restarted by mounting the recorded tape on unit # 8 , specifying 
ITAFE = 2 or 3, and inputting the identical data except for IC which 
must be zero. The following two cards are required for the NFS IBM 

360/67: 

//G0.FTO8FOO1 DD DSN=W0NLIW,UWIT=24OO,V^L=SEE=NFS1OU, 

// DCB=CEECFM=VS ,IiRECL=3204,BLKSIZE=3208 ,DISP= (HEW, KEEP) ,LABEL=( ,SL) 

The boundary conditions are read in on cards. If the shell does 
not have a pole at the first station, XBGINL ^ 0, and cards 7- 15 
describe the boundary conditions at the first station. However, if 
the shell does have a pole at the first station, IBCINL < 0, and 
cards 7-15 are omitted . Cards 7.15 have the format 4El6.8~a.nd correspond 
to the boundary conditions as follows : 


Card 7,16 

Card 8,17 

n(i,i) 

0(1,2) 

0(2,1) 

0(2,2) 

5(3,1) 

5(3,2) 

0(4,1) 

5(4,2) 


Card 9,18 

Card 10,19 

5(1,3) 

0(1,4) 

5(2,3) 

0(2,4) 

5(3,3) 

5(3,4) 

5(4,3) 

5(4,4) 
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Card 11,20 

Card 12,21 

Card 13,22 

Card 14,23 


Card 15,24 

A(l,l) 

A(l,2) 

A(l,3) 

A(l,4) 


u 


4(1)“ 


A(2,l) 

A(2,2) 

A(2,3) 

A(2,4) 


V 


4(2) 


A(3,l) 

A(3,2) 

A(3j3) 

A(3,4) 


w 

— 

i(3) 


A(4,l) 

A(4,2) 

A(4,3) 

A(4,4) 


M 

<3 


4(4) 



If the shell does not have a pole at the final station, IBCFNL £ 0, 
and cards 16-24 describe the boundary conditions at the final station. 
The format and correspondence are the same as for the boundary 
conditions at the first station given above. However, if the shell 
does have a pole at the final station, IBCFNL < 0, and cards 16-24 
are omitted . Note that the boundary conditions are on the total 
variables and not on the individual modes. Thus, it is not possible 
to have different boundary conditions for each mode -without modifying 
the program. An example of the modifications required to change i 
is given in reference 5. Furthermore, note that the boundary conditions 
are input in dimensional form. 

User -Prepared Subroutines 


The geometry of the shell, the inplane and bending stiffnesses 
of the shell, the applied pressure and thermal loads, and the initial 
conditions are introduced to the program through the use of the five 
subroutines GE0M, BDB (K, B, DB, D, DD), FL0 AD(k), TL0AD(K) and INITL. 
This section describes each of these subroutines . 


1. GE0M 

The nondimens ional quantities A, p, y, Ug, u> g , du) s /d§, and p, are 

defined in GE0M as a function of the meridional station number K. 

The correspondence between the nondimens ional variables and the 
F0KTRAN variables is as follows : 


MASS(K) = (n) K = 


h E T 
o o o 


ds J K 

(JWS). 


DEL = A = (meridian length) /[a (KMAX - l)] 
R(K) = (p) K = (r/a) K 

«M(K, . ( V ) K - (ifS) K - (^) K 
0MT(K) = (u> 0 ) K = (a/R 0 ) K 

0MXI(K) = (u) s ) K = (a/R g ) K 

- ( J> K ■ c 2 


y K = 1,2, ... KMAX 


34 



The statements for the static example are 


del=(io5./6.)/iooo. 

D0 4 K=2 ,KMAX 
RK = K 

THET = (EK-!. )*DEL 
R(K) = SIN(THET) 

GAM(K) = C0S (THET)/r(K) 

0MT(K)=1. 

0MXl(K)-l. 

4 DE0M£(K) = 0. 

R(l) = 0. 

GAM(l) = 0. 

0MT(l)=l. 

0MXl(l)=l. 

DE0MX(l) = 0. 

The statements for the dynamic example are 
AKX=KMAX-1 

del=i./akx 

THET=ARSIN(2 .2801/CHAR) 

D 0 11 K=1,KMAX 
AK=K 

R ( K) = ( 7 . 9499+ ( AK-1 . ) * (2 . 2801 ) /AKX ) /CHAR 
GAM(K)=(2 .280l/CHAR) / r(K) 

0MXZ(K)=O. 

DE0MX(K)=O. 

0MT ( K) =€0S ( THET ) /R ( K) 

mass(k)=i. 

11 C0NTXNUE 

2. BDB (K, B, DB, D, DP) 

The nondimens ional stiffness quantities b, db/d|, d, and dd/dc 
are defined in BDB for each meridional station. The correspondence 
"between the stiffness quantities at the Kth station and the FORTRAN 
variables is as follows: 
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B 

DB 

D 

m 


(b) K - (B) K /(E 0 h o ) 

/db% _ /dB^ / a 

l d^K “ ds K E h 
o o 

(d) K = (D) K /(E 0 h3) 

/dd.N _ /dDd / a 
WK _ ds K E h 3 


•) 


) 


where 

B = jEdC/(l-v 2 ) 

D = J*C 2 EdC/(l-v 2 ) 

The statements for the static example are 


b=27.3E+o6/(30,e+o6*i.*(i.-.09)) 

D=B/l2. 

DB=0 . 

DB=0. 


The statements for the dynamic example are 

B=l. 089082 
D=. 09075683 
DB=0. 


DD=0. 


3 . PL0AD(K) 


The nondimens ional Fourier coefficients of the meridional,/ > 
circumferential and normal components of the pressure load, p ^' 11 9 

and - respectively, are defined in PL0AD for each meridional 

station as a function of the Fourier index. In addition, the array 
of Fourier integer numbers n is defined here. The relationship 
between these quantities at the Kth station and F0RTBAN variables is 
as follows: 


HN(M) a n 

rx(M) - (P< n \ - <l s (n) ) K (a/o o h 0 ) 

IT( M ) = (p^\ - (l 9 (n) ) K (a/d 0 b 0 ) 
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M = 1, 2. 


MNMAX 


Hl(M) = (P (n) ) K = (<l (n) ) K ( a /tf 0 h Q ) 


Note that these are stored as inunctions of M only. 

The statements for the static example are 

NN(l)=0 

m(2)=i 

m(3)=3 

PR(1)=-15. 

PR(2)=-19.1 

ER(3)=6.37 

No statements are required for the dynamic example , 
of mode numbers is included in the subroutine INITL. 

4. TL0AD(K) 


The array 


The nondimens ional Fourier coefficients of the thermal loads 

t T^ 5 > if and if ( m T ^ n ' ) ) are defined in TL0AD(K) 

for each meridional station as a function of the Fourier index. The 
FORTRAN variables are defined as follows : 


(n^ 


TT(M) = (t^ VI1 0 K - (e T (n) ) K /(a n h n ) 
d fJ _ (n) u _ / de£ 


o o 


DT(M) = (t T 


)) = (- 
' ; K K ds 


Cn). 




EMT(M) = (m T (n) ) K = ( H T (n) ) K ^ a / CT o h o 3) 

di 

DMT(M) - (^(« T (n) )) K - (-is-) (a 2 /a 0 h o 3 ) 


Note that these are stored as functions of M only. If only thermal 
loads are .applied the array of Fourier inter ger numbers can be 
introduced in TL0AD(K) instead of PL0AD(K) . 


5- INITL 


This subroutine introduces the initial conditions of the non- 
dimensional solution vector z for all the stations, including the 
ficticious stations off the ends of the shell, and all the modes. 
The F0RTRAN variables are defined as follows: 
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Z(I,L) « (z (n) ) K 


U^ n ^(E /acj ) 
v o { cr 

V M (E o /aa 0 ) 

W (n) (E>To> 

_M< n) (a/„ o h o 3)_ 


ZD0T(l,L) 


/dz 


(n) 


dt 


'K 


— r p ^ 

“ o dT 


U (n) (E 0 /M 0 ) 

V (n) (E 0 /aa 0 ) 

W (n) (E 0 /ao 0 ) 

(a/ a h ~ > ') 
s K ’ o o ' 




I 

L 


K 




1 , 2 , 3,4 

1 , 2 , ... 

(KMAX+2)*(MNMAX) 


The index L rims from 1 to KMAX+2 for NW(l) , and from l+KMAX+2 to 
9(KMAX+2) for NN(2), etc. The first element for each value of n 
corresponds to the initial ficticious station, the next element 
corresponds to the first station on the shell, etc. 

The statements for the dynamic example are 

M(l)=0 

MN(2)=1 

1 OT( 3)=2 

NN(4)=4 

PI=3.l4l59 

D0 2 M=1,MAXM 
IT(M.EQ.l) VKL=-444.08/PI 
xf(m.eq.2) vel=-444.o8/2. 

IF(M.EQ,.3) vel=-444 . 08*2 . / ( 3 . *PI ) 
if(m.eq.4) vel=444 . 08*2 . / (15 . *pi ) 

D0 2 K=2,KL 
I=K+l+(M-l)*KMAX2 

2 ZD^T(3,I)=VEL*EMST*TEE0/(CHAR*SIG0)*1O. 


in which KL = KMAX-1 and KMAX2 = KMAX+2. 
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Output Format 


__ _ The output from the program consists of the boundary conditions 

Q, A, and Z at each end of the shell; the input parameters, such as 
the number of stations, the number of modes, etc,; and the cir- 
cumferential coordinates where a summed solution is desired. The 
remainder of the output can appear in either dimensional or non- 
dimensional form. The correspondence between the printed F0RTRAN 
variables and the dimensional and nondimens ional dependent variables 
is given below. For the shell geometry, the following are printed 
at each station: 


RADIUS 

r 

or 

p 

GAMMA 

dr/ds 

r 

or 

dp/p§ 

P 

0MEGA S 


or 

0) 

s 

0MEGA THETA - 


or 

“0 

DE0MEGA S 

& (i /v 

dm 

s 

or -d? 

MASS 

JmdC 

or 

M* 


For the inplane and bending stiffnesses, the following are printed 


at each station: 




B STIFFNESS 

- B 

or 

b 

D STIFFNESS 

- D 

or 

d 

B PRIME 

dB 

ds 

or 

db 

d? 

D PRIME 

cLD 

ds 

or 

dd 

d§ 


For the pressure and thermal loads, the following are printed 
at each station for each value of n for the static analysis: 

N - Fourier index n 


PR 

- 

9 w 

or 

pW 

PK 

- 

^•s 

or 


PT 

- 

>) 

or 

Pa' 1 ” 

TT 

- 

e (n) 

e T 

or 

+ (n) 
t T 
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- »i n) 


(n) 

m T 

MT 

or 


a4 n) 


atM 

DTT 


or 


ds 

dS 


dH^ n) 


dml n) 

DMT 


or 


ds 

d§ 


The following initial conditions are printed at each station for 
each value of n for non zero initial conditions in a dynamic analysis: 


u 


„(a) 

( n ) 

or u v 

V 


v< n > 

(n) 

or v' 

w 

— • 

w (n) 

(n) 

or w v 

M 

s - 


(n) 

or nr 

U 

D)6t 

s 

s 

dl/ n VdT 
dV^ n ^/ dT 
dW^ n VdT 
dM^ n ^/ dT 

V 

d/>t 

- 

W 

dj6t 

- 

M 

s d/3t 

- 


or aJ n Vat 
or dv^ n )/ dt 
or dw^^/dt 
or dm^ n ^/ dt 


Note that station 1 refers to the station on the initial boundary and 
not the ficticious station. 


Every IPRINTth solution is printed with the corresponding load 
factor and the number of iterations . The definitions of the 
printed quantities preceeding the solution are: 


L0AD STEP NUMBER 
TIME STEP NUMBER 
L^AD FACT0R 

TIME 

ITERATIONS 


the number of load intensities considered. 

the number of time steps taken. 

the proportion of the loads given in PL0AD, 

TL0AD and i currently on the shell. 

both nondimens ional and dimensional time are 
given. 

the number of iterations required for convergence. 


The correspondence between the printed terms and the dimensional 
and nondimens ional forces, moments, displacements and rotations is 
as follows : 


N S 

N THETA 
N STHETA 


N or 
s 

N 0 or 
N s0 ° r 



4o 



Q S 


% 

or 

r 

s 

M S 

- 

M s 

or 

m 

s 

M THETA 

- 

g 

CD 

or 

m 0 

M STHETA 

- 

M S0 

or 

m s0 

U 

- 

u 

or 

u 

V 


V 

or 

V 

W 


w 

or 

w 

PHI S 



or 


PHI THETA 

- 

CD 

or 

*0 

PHI 

- 

§ 

or 

tp 


Sample Solutions 

The printed input data and solution for the static analysis 
example is given in figure 6 for the load factor . 744 . The load- 
displacement plot is given in figure 7 for the displacement at the 
pole (station l) and station 3 for 0=0°. The printed input data 
and solution for the dynamic analysis example is given in figure 8 
for T = 500 psec. The time history of the normal displacement at 
s = 6.5 in. and 0 =0° is given in figure 9. 
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— PROBLEM NUMBER 10- 


SAMPLE PROBLEM 10 - IJNSYMMETRICALLY LOADED SPHERICAL CAP, TEST CASE. 

--INPUT DATA RECORD— 

THE BOUNDARY CONDITIONS ARE: 

THE SHELL HAS AN INITIAL POLE 
AT THE FINAL EOGE 

OMEGA BAR LAMBDA BAR EL 


O.L 

0.0 

n.o 

0.0 ) 

N S 

( G.100E Ul 0.0 

0.0 

0.0 

) u 

0.0 

0.0 

u.o 

O.G ) 

N ST + 

( e.o 

0.100E 

01 0.0 

0.0 

) V 

c,<: 

0.0 

0.0 

0.0 ) 

0 S 

( c.o 

0.0 

0. 1 OOF 01 0.0 

) w 

0. 0 

0.0 

0.0 

u.iooe on 

PHI S 

( 0.0 

0.0 

0.0 

0.0 

) M S 


NUMBER OF STATIONS 7 

10 NUMBER OF MODES 4 

INCREMENTAL LOAD FACTOR O.2C0 

MAXIMUM NUMBER OF LOAD STEPS 99 

MAXIMUM NUMBER OF ITERATIONS BQ 

MAXIMUM NUMBER OF LOAD FACTOR CHANGES 2 

CONVERGENCE CRITERION O.ulG 


CHARACTERISTIC SHELL DIMENSION '.MU-JE 

REFERENCE THICKNESS C.lOhlE 01 

REFERENCE ELASTICITY D.30G0E 08 

REFERENCE STRESS 1C; DOE 04 

POISSON’S RATIO 0. 300 JE OC 


CIRCUMFERENTIAL COORDINATES FOR PRINT RECORD, IN RADIANS, ARE: 

0.0 , 


THE DATA IS IN DIMENSIONAL FORM 


Figure 6. Input Data and Solution for the Static Analysis Example 



STATION 

RADIUS 

GAMMA 

OMEGA S 

OMEGA THETA 

DEOMEGA S 

1 

0.0 

c-.o 

0. 1U00E-02 

0 • 1000 E-0 2 

0.0 

2 

0.17 506 32 

C.571AF-°1 

0. 1 00UE-02 

0. 1000E-C2 

0.0 

3 

0.3A99F 02 

C.2856E-01 

0.1000E-02 

0.1000E-02 

0.0 

A 

0. 52 ABF 02 

0. 1903F-01 

0.10G0F-C2 

0.1UGOE-02 

0.0 

5 

C.6994E ^2 

0. 1A26E-01 

O.l'iOOF -02 

G.1000E-C2 

0.0 

6 

0. 0739E 02 

U.11ACF-01 

C.100CE-02 

0.1000E-02 

0.0 

7 

U. JO ABE 03 

U.9A89F-02 

0. 1009E-9? 

C.lOOuE-02 

0.0 


STATION 

B STIFFNESS 

D STIFFNESS 

B PRIME 

0 PRIME 

1 

0.300000E O p - 

0.25000'JE C7 

0.0 

0.0 

2 

n . 300 CODE r 8 

0.250000E 07 

0.0 

0.0 

3 

0.300 00 OE 08 

0. 250000E 07 

0.0 

0.0 

A 

C.3D0C00E 08 

0.250000E 07 

0.0 

0.0 

5 

0. 300C0OE 00 

0. 250000E 07 

0.0 

0.0 

6 

0.300000E 18 

0.250000E 07 

0.0 

0.0 

7 

0. 3 Of) COL* E OP 

C.250000E 07 

o.c 

0.0 



PRESSURE 

AND TEMPERATURE 

COEFFICIENTS 

FOR N= 0 

FOLLOW 



STATION 

PR 


PX 

PT 

TT 

MT 

OTT 

DMT 

1 

-J.150GE 

02 

0.0 

O.C 

O.C 

0.0 

0.0 

0.0 

2 

-0.1500E 

02 

■1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3 

-0. 15GUE 

0? 

u.O 

0.0 

0.0 

0.0 

0.0 

0.0 

A 

-U.1500E 

02 

0.0 

U.O 

0.0 

0.0 

0.0 

0.0 

5 

-0. 1500E 

02 

0. 0 

0.0 

0.0 

0.0 

0.0 

0.0 

6 

-0. 150GE 

02 

o.c 

0.0 

0.0 

0.0 

0.0 

0.0 

7 

-0.1500E 

02 

0.0 

0.0 

o.c 

0.0 

0.0 

0.0 



PRESSURE 

AND TEMPFRATURE 

COEFFICIENTS 

FOR N= 1 

FOLLOW 



STATION 

PR 


PX 

PT 

TT 

MT 

DTT 

DMT 

1 

-C.191GE 

02 

C.L 

0.0 

0.0 

0.0 

0.0 

0.0 

2 

-0. 1910E 

02 

O.c 

0.0 

0.0 

C.O 

0.0 

o.o 

3 

-0.1910E 

02 

c.o 

0.0 

0.0 

0.0 

0.0 

0.0 

A 

-l'. 191 OE 

02 

o.r 

o.c 

0.0 

0.0 

0.0 

0.0 

5 

-0. 191 OE 

02 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

6 

-C.1910E 

02 

'I.C 

0.0 

0.0 

0.0 

0.0 

0.0 

7 

-C. 191 OE 

02 

o.c 

0.0 

c.o 

0.0 

0.0 

0.0 



PRESSURE 

AND TEMPERATURE 

COEFFICIENTS 

FOR N= 3 

FOLLOW 



STATION 

PR 


PX 

PT 

TT 

MT 

DTT 

DMT 

1 

U.637GE 

01 

0.0 

U.O 

0.0 

0.0 

0.0 

0.0 

2 

0.*37PE 

Ul 

c.o 

u.O 

c.o 

0.0 

0.0 

0.0 

3 

C.A37CE 

01 

0.0 

o.c 

0.0 

0.0 

0.0 

0.0 

A 

C.637CE 

01 

U.O 

O.c 

0.0 

0.0 

0.0 

0.0 

5 

C. 637 OE 

01 

o.o 

0.0 

0.0 

0.0 

0.0 

0.0 

6 

0.6370E 

01 

c.o 

0.0 

0.0 

0.0 

0.0 

0.0 

7 

' . 6370E 

01 

■'.o 

0.0 

0.0 

0.0 

0. J 

0.0 


Figure 6. Continued 


THE LOAD STEP NUMBER IS 9 


THE LOAD FACTOR IS C.7440E 00 


THE SOLUTION CONVERGED IN 10 ITERATIONS 


THE SUMMED FORCES, MOMENTS, DISPLACEMENTS AND ROTATIONS FOLLOW FOR THETA = 0.0 


STATION 

N S 

N THETA 

N STHETA 

0 S 

M S 

M THETA 

M STHETA 

1 

- 0. 7200E 04 

-0.4576E 04 

U. 0 

-U.1038E 03 

0.2211E 03 

0.1303E 03 

0.0 

2 

-0.1384F 05 

-0.1361E 0 5 

0.0 

-0. 61 81 E 02 

-0.1626E 04 

-0.1565E 04 

0.0 

3 

-0.1517E 05 

-D.2C11E 05 

0.0 

-0.3293E 02 

-U.2913E 04 

-0.1930E 04 

0.0 

A 

-0.1369E C5 

-0.1995E 05 

0.0 

0.5189E 02 

-0.261 3E 04 

-0.1828E 04 

0.0 

5 

-0.9044E 04 

-0.1389E 05 

0.0 

0.1127E 03 

-0.1166E 04 

-C.1205E 04 

0.0 

6 

-0.5843F 04 

-C.5R93E 04 

0.0 

0.1765E 03 

0.9917E 03 

-0.2303E 03 

0.0 

7 

-0.1257E 04 

-0.3771E 03 

0.0 

0.2868E 03 

0.4425E 04 

0.1328E 04 

0.0 

STATION 

U 

V 

w 

PHI S 

PHI THETA 

PHI 


1 

-0.2584E-01 

0.0 

-0.2880F 00 

0.2900E-01 

0.0 

0.0 


2 

-0.2947E-01 

0.0 

-0.7308E 00 

0.2082E-01 

0.0 

0.0 


3 

-0.2054E-01 

0.0 

-0. 1018 E 01 

0.7398E-02 

0.0 

0.0 


4 

-0.6615F-C2 

0.0 

-0.99O4E 00 

-0.9517E-02 

0.0 

0.0 


5 

0.2747E-02 

o.c 

-U.6R49E 00 

-0.2055E-01 

0.0 

0.0 


6 

0.3677E-G2 

c.e 

-0.2711E 00 

-0.1957E-01 

o.c 

0.0 


7 

-0.5960E-C8 

0.0 

-C.3517E-06 

-0.5058F-08 

0.0 

0.0 






MUDAL OUTPUT FOR 

MODE N = 0 

FOLLOWS 



STATION 

N S 

N THETA 

N STHETA 

Q S 

M S 

M THETA 

M STHETA 

1 

-G.5888E 04 

-0.5888E 04 

0.0 

0.0 

0.1757E 03 

0.1757E 03 

0.0 

2 

-0.6410E 04 

-0.6353E Q4 

0.0 

“0.1433E 02 

0.1757E 03 

0.1222E 03 

0.0 

3 

-0.7064E 04 

-0.6914E 04 

0.0 

-0.2719E C2 

-0.4328E C 3 

-0.1567E 03 

0.0 

4 

-0.6727E 04 

-0.6393E C4 

0.0 

0.1607E 01 

-0 . 49 99 E 03 

-0.3003E 03 

0.0 

5 

-i) • 5 578 E 04 

-0.4661E 04 

0.0 

9.2172E C2 

-0. 2436 E 03 

-0. 2577 E 03 

0.0 

6 

— ?• 46G4E 04 

-0.2505E 04 

0.0 

C.4516E C2 

C.2534F 03 

-0.6132E 02 

0.0 

7 

-0.3185E C4 

-0.9554E 03 

0.0 

0.8682E 02 

0.1211E 04 

0.3634F 03 

0.0 

STATION 

U 

V 

w 

PHI S 

PHI THETA 

PHI 


1 

0.0 

0.0 

-0.2880E 00 

U.O 

0.0 

0.0 


2 

-0.291 8E~P 3 

0.0 

-C.2880E 00 

0 . 534 3E-C3 

0.0 

u.o 


3 

0.1396E-02 

J.G 

-C.3J67E Ou 

-0.4142F-03 

0.0 

o.c 


4 

0.3637E-02 

0.0 

-0.2735E Of 

-0.3473E-02 

o.c 

0. 0 


5 

n . 433 8E-G2 

c.u 

-J.18 50E UC 

-0.5689E-0? 

0.0 

O.n 


6 

0.2750F-02 

fj.0 

-0.7421E-01 

-O.5284E-02 

0.0 

0.0 


7 

0.62J9E-09 

U.o 

-0.6358E-07 

U.1R1 7E-C8 

c.o 

o.c 



Figure 6 


Continued 



MODAL OUTPUT FOR MODE N = 1 FOLLOWS 


STATION 

N S 

N THETA 

N STHETA 

0 S 

M S 

M THETA 

M STHETA 

1 

Cj . C 

0.0 

0 • l 

-0.1038E 03 

0. 0 

0.0 

0.0 

2 

-0.6742E C'A 

-0.748 IF 0 4 

-7.I83RE U4 

-0.5941E 02 

-0.1816E 04 

-0.1264E 04 

0.5525E 03 

3 

-U.7326E 04 

-D.1222E '■*> 

-0.4532E 04 

0.185BE 01 

-0.2G80E 04 

-0.1399E 04 

0.5952E 03 

4 

6906F 04 

-0.1234E 05 

-0.5566E 04 

0 • 38 30 F 02 

-C.1665E 04 

-0.1244E 04 

0.5720E 03 

5 

-0* 4 790E 04 

-9.8735E -4 

-0.5549E 04 

0.6816E J2 

-C.8401E 03 

-0.8649E 03 

0.4709F 03 

6 

-0.2723E U 4 

-C.34D7E 14 

-C.47U1E 04 

0. 1 21 RE 03 

0.4724E 03 

— 0. 2416 E 03 

0.2954E 03 

7 

U.6538E 03 

9.1962E 13 

-0.3387E 04 

0.2270E 03 

0.3019E 04 

0.9D58E 03 

-0.2821 E 00 

STATION 

U 

V 

W 

PHI S 

PHI THETA 

PHI 


1 

-J.2584E-01 

0. 2584E-DI 

u.o 

0.290UE-D1 

-G.290UE-01 

O.u 


2 

-0.2584E-01 

0.2870E-H 

-Q.4113F 0 0 

0.1795E-01 

-0. 2347E-01 

0. 8296 E— 04 


3 

-0.1848E-G1 

0.2593F-C1 

-0.6290 E 00 

0.6036E-02 

-0.1795E-01 

-0. 1298E-04 


4 

-1.8680E-C2 

0.2C35F-01 

-0.6232E 00 

-0.5318F-02 

-U.1185E-01 

-U.7C02E-04 


5 

-0.1416E-C2 

0.1 3 27E-C 1 

- 0 . 4 4 3 2 E 0 0 

-0.1252E-C1 

-0.6323E-0 2 

-0. 1 175E-03 


6 

J. 9771E-C3 

0.6213E-L2 

-0.1849E OU 

-0. 1266F-0 1 

-0. 2 1 10E-02 

-0.1485E-03 


7 

-0.5588E-C8 

0.1192E-1-7 

3497E-Q6 

-0.8180E-08 

-0. 3324E-08 

-0.1613E-03 






MODAL OUTPUT FOR 

MODE N = 3 

FOLLOWS 



STATION 

N S 

N THETA 

N STHETA 

0 S 

M S 

M THETA 

M STHETA 

1 

0.0 

0.0 

C .n 

0.0 

0.0 

0.0 

0.0 

2 

0.8360E 03 

-0.6996E 03 

-C.7687E 03 

0.8236E 01 

-0.3132F 02 

0.2614E 03 

0.1267E 03 

3 

9. 121 IE 04 

-0.4578E 03 

-C.8223E 03 

0.4B83E 01 

C.1137E 03 

0.3885E 03 

0.6502E 02 

4 

0 • 1499 E 04 

0.2354E C 3 

-0. 3732E 03 

0.1461E Cl 

0.2192E 03 

0.3771 F 03 

-O.llllE 02 

5 

0.1669E 04 

0.7190E f 3 

9.1172F C3 

-0.6189E Cl 

U.2923E 03 

0.3210E 03 

-0.8573E 02 

6 

0.1137E 04 

0.5551E 03 

0.411 IF 0 3 

-0. 2778 E 02 

0* 1456E C3 

0.1705E U3 

-0.1129E 03 

7 

9.1267E 13 

0.3796E i:2 

0.3961E 03 

-U.6937E 02 

-0.5344E 03 

-0.1603E 03 

0.3286E-01 

STATION 

U 

V 

W 

PHI S 

PHI THETA 

PHI 


1 

G.O 

0.0 

0.0 

U.O 

0.0 

0.0 


2 

-0 • 3575E-04 

0 ■ 6 6 0 4E - L 4 

0. 5479E-U2 

-C.7353E-03 

0.9394 E-03 

-0.2211 E-05 


3 

C.196SE-C3 

-G.7231E-C 4 

^.2S73E-"1 

-U.1168E-02 

0.2206E-02 

-0.1396E-05 


4 

0.4647E-03 

-0.5498E-03 

0.463RE-O1 

-U.7714E-C3 

r, .2651E-02 

-0.4023F-05 


5 

0.31 53E-03 

-U.91 75E-03 

■J. 5275E-01 

0. 39C1E-C3 

0.2262F-02 

-0.1653 E-05 


6 

-D.3925E-C4 

- 0 . 6 8 U 8 E - u 3 

3 2 74 E~ 01 

0.15U7F-02 

C.1123E-C2 

O.0553E-O5 


7 

-J. 1242E-09 

'J.2111E-CP 

u. 63 SbE-07 

C.2214E-C8 

9 . 18 22F-08 

0.1886E-04 
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Figure 7. Load-Displacement Curves for 

Static Analysis Example Problem 





-PROBLEM NUMBER 108 


LMSC TEST CASE, IMPULSIVELY LOAOED CONE 


— INPUT DATA RECORD— 


THE BOUNDARY CONDITIONS ARE: 

AT THE INITIAL EDGE 
OMEGA BAR 


-LAMBDA BAR- 


£ 


AT THE FINAL EDGE 
OMEGA BAR— 


-LAMBDA BAR- 


NUMBER OF STATIONS 31 

NUMBER OF MODES 4 

INCREMENTAL TIME 0.182E-01 

MAXIMUM NUMBER OF TIME STEPS 750 

MAXIMUM NUMBER OF ITERATIONS 10 

CONVERGENCE CRITERION 0.010 

CHARACTERISTIC SHELL DIMENSION 0.1500E 02 

REFERENCE THICKNESS 0.5430E 00 

REFERENCE ELASTICITY 0.3520E 07 

REFERENCE STRESS 0. 1000E DE- 
REFERENCE TIME 0.1097E-03 

POISSON* S RATIO 0.2860E 00 


—EL 


( 0.0 

0.0 

0.0 

0.0 ) 

N S 

( 0.100E 

01 0.0 

0.0 

0.0 

) U 

0.0 

( 0.0 

0.0 

0.0 

0.0 ) 

N ST + 

( 0.0 

0.1C0E 

01 0.0 

0.0 

) V 

0.0 

< 0.0 

0.0 

0.0 

0.0 ) 

Q S 

( C.O 

0.0 

0.100E 01 0.0 

) w 

0.0 

( 0.0 

0.0 

c.o 

0.100E 01) 

PHI S 

( 0.0 

0.0 

0.0 

0.0 

) M S 

0.0 


EL 


0.0 

0.0 

0.0 

0.0 ) 

N S 

( 0.100E 01 

0.0 

0.0 

0.0 

) U 

0.0 

0.0 

0.0 

0.0 

0.0 ) 

N ST + 

( 0.0 

0.100E 

01 0.0 

0.0 

) V 

0.0 

0.0 

0.0 

0.0 

0.0 ) 

Q S 

( 0.0 

0.0 

0.1C0E 

01 0.0 

) w 

0.0 

0.0 

0.0 

0.0 

0.100E 01) 

PHI S 

( 0.0 

O.G 

0.0 

0.0 

) M S 

0.0 


CIRCUMFERENTIAL COORDINATES FOR PRINT RECORD, IN RADIANS, ARE: 
0.0 , 0.314159E 01 , 


THE DATA IS IN DIMENSIONAL FORM 


Figure 8. Input Data and Solution for the Dynamic Analysis Example 



STATION 

RADIUS 


GAMMA 

OMEGA S 

OMEGA THETA 

DEOMEGA S 

MASS 

1 

G.7950E 

01 

0 • 1912E-01 

0.0 

0.1243E 

00 

0.0 

0.1021E-03 

2 

0.8026E 

01 

G.1893E-U1 

0.0 

0.1231E 

00 

0.0 

G.1021E-03 

3 

0# 8102E 

01 

C • 1 876E-01 

0.0 

G.1220E 

00 

0.0 

0.1021E-03 

4 

0* 8178E 

01 

0.1 858F-01 

C'.O 

0.1209E 

00 

0.0 

0.1021E-03 

5 

0.8254E 

01 

0.1841E-01 

0.0 

0.1197E 

00 

O.G 

0. 1021E-03 

6 

0.8330E 

01 

0.1824E-01 

0.0 

0.1187E 

00 

0.0 

0.1021E-03 

7 

0.8406E 

01 

0.1808E-C1 

0.0 

C.1176E 

OG 

0.0 

0.1021E-03 

8 

0.8482E 

01 

0. 1 792E-01 

0.0 

C.1165E 

00 

0.0 

G.1021E-93 

9 

0.8558E 

01 

0.1776E-01 

0.0 

G.1155E 

00 

0.0 

0.1021E-03 

10 

0.8634E 

01 

0.1760E-01 

0.0 

0.1145E 

00 

0.0 

0. 1021E-03 

11 

0.8710E 

01 

0.1745E-01 

0.0 

0.1135E 

00 

0.0 

0. 1021E-03 

12 

0.8786E 

01 

0.1 730 E -01 

0.0 

0.1125E 

00 

0.0 

0.1021E-03 

13 

0.8862E 

01 

0.1715E-01 

0.0 

G.1115E 

oc 

0.0 

0. 1021E-03 

14 

0.8938E 

01 

0.1700E-01 

0.0 

0.1106E 

00 

0.0 

0.1021E-03 

15 

0.9014E 

01 

0.1686E-01 

0.0 

0.1097E 

00 

0.0 

C.1021E-03 

16 

0.9090E 

01 

0.1 672E-01 

0.0 

0.1087E 

00 

0.0 

G.1021E-03 

17 

0.9166E 

01 

0.1658E-01 

0.0 

0.1078E 

00 

0.0 

0.1021E-03 

18 

0.9242E 

01 

0.1644E-01 

0.0 

0.1069E 

00 

0.0 

0. 1021E-03 

19 

0.9318E 

01 

0.1631E-01 

0.0 

0.1061E 

00 

0.0 

0.1021E-03 

20 

0.9394E 

01 

C.1618E-01 

0.0 

0.1052E 

00 

0.0 

0.1021E-03 

21 

0.9470E 

01 

0. 1 605E-01 

0.0 

0.1044E 

00 

0.0 

0. 1021E-03 

22 

0.9546E 

01 

0. 1 592E-01 

0.0 

0.1035E 

00 

0.0 

0.1021E-03 

23 

0.9622E 

01 

0*1 579E-01 

0.0 

0.1027E 

00 

0.0 

C.1021E-03 

24 

0.9698E 

01 

0.1567E-01 

0.0 

0.1019E 

00 

0.0 

0.1021E-03 

25 

0.9774E 

01 

0.1 555E-01 

0.0 

0. 10 HE 

00 

0.0 

0.1021E-03 

26 

0.9850E 

01 

0.1543E-01 

0.0 

0.1003E 

oc 

0.0 

0.1021E-03 

27 

0.9926E 

01 

0.1531E-01 

0.0 

0.9958E-01 

0.0 

0. 1021E-03 

28 

0.1000E 

02 

0.1519E-01 

0.0 

0.9882E- 

-01 

0.0 

0. 1021E-03 

29 

0.1C08E 

02 

0. 1508E-01 

0.0 

0.9807E-0 1 

0.0 

0. 1021E-03 

30 

0.1015E 

02 

0.1497E-01 

0.0 

0.9734E-01 

0.0 

0. 1021E-03 

31 

0.1023E 

02 

0.1485E-01 

0.0 

0.9662E- 

-01 

0.0 

0.1021E-03 


ST AT ION 

B STIFFNESS 

D STIFFNESS 

B PRIME 

D PRIME 

1 

0.208163E 

07 

0. 51 1471 E 

05 

0.0 

0.0 

2 

0.208163E 

07 

0.511471E 

05 

o.c 

0.0 

3 

0. 2081 63 E 

07 

0.511471E 

05 

0.0 

0.0 

4 

0.208163E 

07 

0. 511471 E 

05 

0.0 

0.0 

5 

0.2C8163E 

07 

0. 511471E 

05 

O.G 

0.0 

6 

0.208163E 

07 

0.51 1471 E 

05 

0.0 

0.0 

7 

0.208 1 63E 

07 

0. 511471 E 

05 

0.0 

0.0 

8 

C.298163E 

07 

0.51 1471 E 

C5 

o.c 

0.0 

9 

0.2G8163E 

07 

0. 511471 E 

05 

0.0 

c.o 

1C 

0.208163E 

07 

0.51 1471E 

05 

o.c 

0.0 

11 

0.208163E 

07 

0.51 1471 E 

05 

o.c 

O.G 

12 

0.2C8163E 

07 

0.511471E 

C 5 

o.c 

0.0 

13 

0.2C8163E 

07 

C.511471E 

05 

O.G 

o.c 

14 

0.2C 8163E 

07 

0.51 1471 E 

05 

{1.0 

0.0 

15 

0.208163E 

07 

0.511471E 

C5 

0.0 

0. J 

16 

0.208163E 

07 

0.51 1471 E 

C 5 

o.c 

0.0 

17 

0.2G8163E 

07 

0. 51 1471 E 

05 

0.0 

o.c 

18 

0. 2°8 1 63E 

07 

0.51 1471 E 

u5 

0.0 

0.0 

19 

0.208163E 

07 

0. 511471 E 

05 

o.c 

0.0 

20 

0.208163E 

07 

0. 51 1471 E 

05 

0.0 

O.G 

21 

0.208163E 

07 

U.511471E 

05 

O.f; 

0.0 

22 

0. 2081 63 E 

07 

0. 511471 E 

05 

0.0 

0.0 

23 

0.208163E 

07 

0.51 1471 E 

05 

0.0 

0.0 

24 

0.208 1 63 E 

07 

0. 51 1471 E 

05 

o.c 

0.0 

25 

0.2C81 63E 

C7 

0. 511471 £ 

C5 

O.G 

0.0 

26 

0. 208 i 63 E 

07 

0.511471E 

05 

O.C; 

0.0 

27 

0.208163E 

07 

U. 511471 E 

05 

C . 0 

c.o 

28 

0.208163E 

07 

0. 511471 E 

C 5 

0.0 

0.0 

29 

0.2G8163E 

07 

G.511471E 

05 

0. o 

0.0 

3u 

0.2G8163E 

07 

0. 51 1471 E 

05 

C.u 

o.c 

31 

0.2C-8163E 

07 

0. 51 1471 E 

05 

c.o 

0.0 


Figure 8. Continued 



THE INITIAL CONDITIONS FOR N= 0 FOLLOW 


STATION 

U 

V 

w 

M S 

1 

0.0 

0.0 

0.0 

0.0 

2 

o.c 

o.c 

o.c 

0.0 

3 

0.0 

0.0 

0.0 

0.0 

4 

0.0 

0.0 

0.0 

0.0 

5 

0.0 

0.0 

0.0 

0.0 

6 

0.0 

0.0 

0.0 

0.0 

7 

0.0 

0.0 

o.c 

0.0 

8 

0.0 

0.0 

0.0 

0.0 

9 

0.0 

0 . c 

0.0 

0.0 

10 

0.0 

0.0 

0.0 

0.0 

11 

0.0 

0.0 

o.c 

0.0 

12 

0.0 

0.0 

0.0 

0.0 

13 

o.c 

0.0 

0.0 

0.0 

14 

0.0 

0.0 

o.c 

0.0 

15 

0.0 

0.0 

0.0 

0.0 

16 

0.0 

0.0 

0.0 

0.0 

17 

o.c 

o.c 

0.0 

o.c 

18 

o.c 

0.0 

0.0 

0.0 

19 

0.0 

0.0 

o.c 

0.0 

20 

o.c 

0.0 

0.0 

0.0 

21 

0.0 

0.0 

o.c 

0.0 

22 

0.0 

0.0 

o.c 

0.0 

23 

0.0 

o.c 

G.O 

0.0 

24 

O.L 

0.0 

0.0 

0.0 

25 

0.0 

0.0 

0.0 

0.0 

26 

o.c 

o.c 

0.0 

0.0 

27 

o.c 

0.0 

0.0 

0.0 

28 

o.c 

0.0 

0.0 

0.0 

29 

0.0 

0.0 

0.0 

0.0 

30 

0.0 

o.c 

o.c 

0.0 

31 

0.0 

0.0 

o.c 

0.0 


STATION 

U DOT 

V DOT 

W DOT 


M S DOT 

1 

O.G 

0.0 

c.o 


0.0 

2 

0.0 

0.0 

-0. 141355E 

04 

0.0 

3 

0.0 

o.c 

-0.141355E 

04 

0.0 

4 

0.0 

0.0 

-0.141355E 

04 

0.0 

5 

0.0 

0.0 

-0. 141355E 

04 

0.0 

6 

O.C 

0.0 

-0.141355E 

04 

0.0 

7 

O.C 

0.0 

-0. 141355 E 

04 

0.0 

8 

0.0 

0.0 

-0. 141355E 

04 

0.0 

9 

0.0 

0.0 

-0. 141355E 

04 

0.0 

1C 

0.0 

0.0 

-0.141355E 

04 

0.0 

11 

O.Q 

0.0 

-0.141355E 

04 

0.0 

12 

0.0 

0.0 

-0. 141355E 

04 

0.0 

13 

0.0 

0.0 

-G.141355E 

04 

0.0 

14 

0.0 

0.0 

-0.141355E 

04 

0.0 

15 

o.c 

0.0 

-0. 141355E 

04 

0.0 

16 

0.0 

0.0 

-0. 141355E 

04 

0.0 

17 

0.0 

0.0 

-0. 141355E 

04 

0.0 

18 

0.0 

0.0 

-0.141355E 

04 

0.0 

19 

0.0 

0.0 

-0.141355E 

04 

0.0 

20 

o.c 

0.0 

-0. 141355E 

C4 

0.0 

21 

0.0 

0.0 

-0.141355E 

04 

0.0 

22 

0.0 

0.0 

-0. 141355E 

04 

0.0 

23 

0.0 

0.0 

-0.141355E 

04 

0.0 

24 

0.0 

0.0 

-0. 1413 55 E 

04 

0.0 

25 

o.c 

o.c 

-0. 141355E 

04 

0.0 

26 

0.0 

0.0 

-0.141355E 

C4 

0.0 

27 

0.0 

0.0 

-0. 141355 E 

04 

0.0 

28 

0.0 

0.0 

-0. 141355E 

04 

c.o 

29 

o.c 

o.c 

-0.141355E 

04 

0.0 

30 

O.G 

o.c 

-0. 141355E 

04 

0.0 

31 

0.0 

0.0 

0.0 


0.0 


Figure 8. Continued 


THE INITIAL CONDITIONS FOR N= 1 FOLLOW 
STATION U 


I 

O.C 

2 

0.0 

3 

0.0 

4 

0.0 

5 

O.C 

6 

0.0 

7 

0.0 

8 

0.0 

9 

0.0 

IC 

0.0 

11 

0.0 

12 

0.0 

13 

0.0 

14 

0.0 

15 

0.0 

16 

0.0 

17 

0.0 

18 

0.0 

19 

0.0 

20 

0.0 

21 

0.0 

22 

0.0 

23 

0.0 

24 

0.0 

25 

0.0 

26 

0.0 

27 

0.0 

28 

0.0 

29 

0.0 

30 

0.0 

31 

0.0 


STATION 

U DOT 

1 

0.0 

2 

0.0 

3 

0.0 

4 

0.0 

5 

0.0 

6 

0.0 

7 

O.C 

8 

0.0 

9 

0.0 

10 

0.0 

11 

0.0 

12 

0.0 

13 

0.0 

14 

O.C 

15 

0.0 

16 

0.0 

17 

O.C 

18 

0.0 

19 

0.0 

20 

O.C 

21 

0.0 

22 

0.0 

23 

O.C 

24 

0.0 

25 

0.0 

26 

0.0 

27 

0.0 

28 

0.0 

29 

0.0 

30 

O.C 

31 

O.C 


Figure 8. Continued 


V 

w 

M s 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

c.c 

O.G 

0.0 

0.0 

0.0 

0.0 

U.O 

O.G 

0.0 

0.0 

0.0 

O.G 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

c.o 

0.0 

0.0 

0.0 

0.0 

0.0 


V DOT 

W DOT 


M S DOT 

0.0 

0.0 


0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

“ 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0. 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0.2220 40 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

O.C 

- 0 . 2220 40 E 

04 

0.0 

O.C 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0. 222040 E 

C 4 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

O.C 

- 0 . 22204 GE 

04 

0.0 

0.0 

- C . 222040 E 

04 

0.0 

0.0 

- 0.222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0. 222040 E 

04 

0.0 

0.0 

- 0 . 222040 E 

04 

0.0 

0.0 

- 0. 222040 E 

04 

0.0 

O.C 

0.0 


0*0 


THE INITIAL CONDITIONS FOR N= 2 FOLLOW 


STATION 

U 

1 

0.0 

2 

0.0 

3 

0.0 

4 

o.c 

5 

o.c 

6 

0.0 

7 

0.0 

8 

0.0 

9 

0.0 

10 

0.0 

11 

0.0 

12 

0.0 

13 

0.0 

14 

0.0 

15 

0.0 

16 

0.0 

17 

0.0 

18 

0.0 

19 

0.0 

20 

0.0 

21 

0.0 

22 

0.0 

23 

0.0 

24 

0.0 

25 

o.c 

26 

0.0 

27 

0.0 

28 

0.0 

29 

0.0 

30 

0.0 

31 

o.c 



V 


W 



V DOT 

W DOT 


M S DOT 

0.0 

0.0 


0.0 

0.0 

-0.942367E 

C 3 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0. 0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0. 942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

o.c 

-C.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0. 942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

C . 0 

-0. 942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

o.c 

-0.942367E 

C3 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

o.c 

-0.942367E 

03 

0.0 

C • 0 

-C.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-C.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

o.c 

-C.942367E 

03 

0.0 

o.c 

-0.942367E 

03 

0.0 

0. G 

-C.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

-0.942367E 

03 

0.0 

0.0 

O.C 


0.0 



THE INITIAL CONDITIONS FOR N = 4 FOLLOW 


STATION 

U 

V 

w 

M S 

1 

0.0 

0.0 

c.o 

0.0 

2 

0.0 

o.c 

0.0 

0.0 

3 

o.c 

0.0 

0.0 

0.0 

4 

0.0 

0.0 

0.0 

0.0 

5 

o.c 

0.0 

0.0 

0.0 

6 

o.c 

0.0 

0.0 

0.0 

7 

0.0 

o.c 

0.0 

0.0 

8 

0.0 

0.0 

0.0 

0.0 

9 

0.0 

0.0 

0.0 

0.0 

10 

0.0 

0.0 

0.0 

0.0 

11 

0.0 

0.0 

0.0 

0.0 

12 

o.c 

0.0 

0.0 

0.0 

13 

0.0 

0.0 

0.0 

0.0 

14 

0.0 

0.0 

0.0 

0.0 

15 

0.0 

0.0 

0.0 

0.0 

16 

0.0 

o.c 

0.0 

0.0 

17 

0.0 

0.0 

0.0 

0.0 

18 

o.c 

0.0 

0.0 

0.0 

19 

0.0 

0.0 

0.0 

0.0 

20 

0.0 

o.c 

0.0 

0.0 

21 

0.0 

0.0 

0.0 

0.0 

22 

o.o 

c.o 

0.0 

0.0 

23 

0.0 

0.0 

0.0 

0.0 

24 

o.c 

0.0 

o.c 

0.0 

25 

0.0 

0.0 

0.0 

0.0 

26 

0.0 

o.c 

0.0 

0.0 

27 

o.c 

0.0 

0.0 

0.0 

28 

0.0 

0.0 

0.0 

0.0 

29 

0.0 

0.0 

o.c 

c.o 

3C 

o.c 

0.0 

0.0 

0.0 

31 

o.c 

0.0 

0.0 

0.0 


STATION 

U DOT 

V DOT 

W DOT 


M S DOT 

1 

0.0 

0.0 

0. 0 


0.0 

2 

0.0 

o.c 

C.188473E 

03 

0.0 

3 

0.0 

c.o 

C.188473E 

03 

0.0 

4 

0.0 

o.c 

0. 188473E 

03 

0.0 

5 

0.0 

0.0 

0.188473E 

03 

0.0 

6 

0.0 

0.0 

0.188473E 

03 

0.0 

7 

0.0 

o.c 

0. 188473E 

03 

0.0 

8 

0.0 

0.0 

0. 188473E 

03 

0.0 

9 

o.o 

O.D 

0. 188473E 

03 

0.0 

10 

0.0 

0.0 

0. 188473E 

03 

0.0 

11 

o.o 

0.0 

0. 188473E 

03 

0.0 

12 

o.c 

0.0 

0. 188473E 

03 

0.0 

13 

o.c 

0.0 

0.188473E 

03 

0.0 

14 

o.c 

0.0 

0. 188473E 

03 

0.0 

15 

0.0 

0.0 

0.188473E 

03 

0.0 

16 

o.c 

0.0 

0. 188473E 

G3 

0.0 

17 

o.o 

0.0 

0. 188473E 

03 

0.0 

18 

o.c 

0.0 

0. 188473E 

03 

0.0 

19 

o.o 

o.c 

C.188473E 

03 

0.0 

20 

0.0 

0.0 

0. 188473E 

03 

0.0 

21 

o.c 

o.c 

0.1B8473E 

03 

0.0 

22 

0.0 

0.0 

0. 188473E 

03 

0.0 

23 

o.c 

0.0 

0.188473E 

03 

0.0 

24 

0.0 

o.c 

0. 188473E 

03 

0.0 

25 

o.c 

c.o 

0.188473E 

03 

0.0 

26 

0.0 

0.0 

0.188473E 

C 3 

0.0 

27 

o.c 

o.c 

0. 188473E 

03 

0.0 

28 

o.c 

0.0 

0. 188473E 

03 

0.0 

29 

0.0 

0.0 

0. 188473E 

03 

0.0 

30 

o.c 

0.0 

0.188473E 

03 

0.0 

31 

0.0 

o.c 

o.c 


0.0 


Figure 8 . Continued 



THE TIME STEP NUMBER IS 250 


THE TIME IS 4.56 OR 9.50CE-03 SECONDS 


THE SOLUTION CONVERGED IN 2 ITERATIONS 


THE SUMMED FORCES, MOMENTS, DISPLACEMENTS AND ROTATIONS FOLLOW FOR THETA = 0.0 


STATION 

N S 

N THETA 

N STHETA 

Q S 

M S 

M THETA 

M STHETA 

1 

-0.3495E 04 

-0.9659E 03 

0.0 

-0.2371E 04 

0.2368E 03 

0.6771E 02 

0.0 

14 

0.7610E 03 

0.2495E 04 

0.0 

0.2499E 04 

0. 18 38 E 04 

0.4449E U3 

0.0 

27 

0.3077E 04 

-0.3228E 04 

0.0 

-0.2860E 04 

-0.6161E 03 

-0.1717E 03 

0.0 

31 

0.2633E 04 

0.7879E 03 

C.O 

0.1325E 04 

-0.9602E 02 

-0.2746E 02 

0.0 

STATION 

U 

V 

W 

PHI S 

PHI THETA 

PHI 


1 

0.0 . 

0.0 

0.2250E-09 

-0.2540E-09 

0.0 

0.0 


14 

-0.6019E-02 

0.0 

0.1312E 00 

-0.3548E-01 

0.0 

0.0 


27 

-0.4080E-02 

0.0 

0.3925E-01 

0. 4481E-01 

0.0 

0.0 


31 

-0.2044E-09 

0.0 

-0.7940E-10 

-0.5080E-09 

0.0 

0.0 



THE SUMMED FORCES 

, MOMENTS, 1 

DISPLACEMENTS AND 

ROTATIONS FOLLOW 

FOR THETA = 

0 • 314159E 01 



STATION 

N S 

N THETA 

N STHETA 

Q S 

M S 

M THETA 

M STHETA 

1 

0.4581E 04 

C.1276E 04 

-0.2199E-01 

-0.1034E 04 

0.1454E 04 

0.4158E 03 

-0.1008E-03 

0.5865E-03 

14 

-0.6556E 04 

-0.3293E 05 

“0. 8604E-02 

-0.2622E 03 

-0.8920E 03 

-0.3827E 03 

27 

-0.3589E 04 

-0.6110E 04 

0.2565E-01 

0.3636E 03 

0.5828E 03 

0.1073E 03 

-0.4031E-03 

31 

-0.1640E 04 

-0.5054E 03 

0. 1780E-01 

0. 5451 E 03 

0.1496E 04 

0.4279E 03 

0.6339E-04 

STATION 

U 

V 

W 

PHI S 

PHI THETA 

PHI 


1 

0.0 

0.0 

-0. 4845E-09 

0.3302E-08 

0.6000E-15 

-0. 1480E-07 


14 

0.5235E-G2 

-0.2114E-06 

-0.2075E 00 

0. 2528E-01 

C. 98 15E-07 

-0.1485E-07 


27 

0.2690E-02 

-0.6559E-07 

-0.4600E-01 

-0.3938E-01 

0.23 16E-07 

0.2331E-07 


31 

0.8138E-1G 

-0.2562E-14 

“0. 4764E-10 

0. 2540E-08 

-0.2530E-15 

0.1198E-07 



Figure 8. 
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Subroutine Descriptions 


MAIN 


This program controls the logical connections between the subroutines. 
The case description, control parameters, physical constants, and 
boundary conditions are both read and printed out in this routine. The 
boundary conditions are nondimens ionali zed and many of the common 
indices and coefficients are determined here. The iteration procedure, 
the load incrementing procedure, and the calculation for the estimate 
of the next solution are all carried out here. The data for re- 
starting the computation is written on tape and read from tape here. 


Subroutine GE0M 

This subroutine computes the nondimens ional geometry functions 
of the shell. 


Subroutine BDB (K,B,DB,D,DD) 

This subroutine computes the nondimens ional inplane and bending 
stiffnesses of the shell. 


Subroutine PL0 AD(k) 

This subroutine computes the nondimens ional Fourier coefficients 
of the loads applied to the shell. 

Subroutine TL0AD(K) 

This subroutine computes the nondimens ional thermal loads. 


Subroutine INITL 

This subroutine computes the initial conditions on z and dz/dt. 


Subroutine PMA.TRX 

This subroutine calls subroutines HJ(K,MN), EFG(K,MN), ABC, and 
PAM)D(K,MN) to set up the P, (P), P, (DEE), and f, (DST), matrices 
given by equations (30). Matrices DL, DG, and DF are set up for the 
calculation of x 1 given by equation ( 31 a), where 


x x = DLj^ + DGg 1 + DFf 1 
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The special P matrix for a shell with an initial pole, given in Ref. 
[l], is also computed here. Matrices ZF1M, ZF2M, ZF3M, and ZF4 m 
are set up for the calculation of given by equation (31b) 

where 


« ZF1Mj^ k + ZF2M5C k + ZFSMX^ + ZF4Mf K 

If the shell has a final pole, the matrices CLO, CL1 and CIj 2 are 
prepared for the calculation of Z given by equation (D-3) in Ref. 
[l], where 

Z K « CLOX K-1 or CL1X k _ 1 or CISXg-1 


depending upon whether n = 0 , 1 or 2. 


Subroutine HJ(K,MN) 

This subroutine computes the elements of the H and JAY matrices 
for both boundaries of the shell. The elements of H and JAY are 
defined in Ref. [l]. 

Subroutine EFG(K,MN) 

This subroutine prepares the elements of the E, F, and G matrices 
for each meridian station K and for each Fourier mode MN. The 
matrices E, F, and G are given in Ref. [l 3 . 


Subroutine ABC 

This subroutine computes the elements of the A, BEE, and C matrices 
defined by equation (25) . 


Subroutine PANDD(K,MN) 

This subroutine computes the elements of the P, (p) , P, (DEE), 
and £,(DST), matrices for each meridian station K and Fourier mode 
MN. These matrices are computed and saved because they do not change 
during either the iteration procedure or the load increment procedure, 
i.e., they are a function of the shell’s initial geometry and stiff- 
ness . 
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Subroutine XAUDZ 


This subroutine computes the x vector using the P, P, and P 

matrices and solves for the z vector for the applied and pseudo 

loads. The subourtines PHXBET(K) and TEAETA(k) are called and 

the previous solution for z, or the estimated value of z, is used 

to calculate the nonlinear Beta and Eta terms. The matrices FFS 

and FLS are the values of f at the initial and final edges of the 

shell. The subroutine F0 RCE(k) is called to calculate the load 

vector g, (GEE), and the x vector at each meridian station. Once 

the x vector is obtained for all meridian stations the solution 

for z given by equation (31b) is obtained, and the solution 
K+l 

for at all the other meridian stations defined by equation (29b) 

is obtained. The solution z at the imaginary station off the 
initial edge of the shell is obtained last. The test for 
convergence of the solution is made as z is computed. The special 
conditions for computing z at either an initial or a final pole are 
also in this routine. 


Subroutine INLP0L 

This subroutine computes, for a shell with an initial pole, 

the nonlinear terms P , P Q , p fl , T] , and T| Q at the pole. The 

s u su ss us 

appropriate equations are given in Ref. [l]. 


Subroutine FNLP0L 

This subroutine computes, for a shell with a final pole, the 

nonlinear terms p , P Q , P Q , T| , and TU at the pole. The 
s u s u s s us 

appropriate equations are given in Ref. [l]. 


Subroutine ft0DES 

In MjZfoES, arrays that define those sets of indices that combine 
to equal each value of n in the problem are determined. M0DES is 
called prior to the first iteration and after every iteration until 
a specified number of Fourier terms is reached. Each Fourier index 
in the problem is subtracted from all other Fourier indices and the 
result is compared with all Fourier indices to see if the new value 
exists in the program. (The same comparison is never made twice.) 

If it does, the locations of the two indices that made the combination 
are stored in two special two-dimensional arrays, IB and JD. One 
argument of each array is the value of the new index and the other 
is the number of combinations of indices that also give this value 
of the index. If there is no index in the program that matches the 
new one, then a new Fourier term has been generated and will be 
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considered, in the next iteration for solution. The variable MAXD 
stores the total number of such combinations for each value of the 
Fourier index. In a similar manner, each index is added to every 
other index and the sum compared with all indices. This result is 
stored in the two two-dimensional arrays, IS and JS, in the same 
manner as was done for the subtraction case. The variable MAXS 
stores the total number of summation combinations for each value 
of the Fourier index. A special routine handles the cases where the 
index is added to and subtracted from itself. The two-dimensional 
array IJS stores the location of the index and the variable MAXSY 
stores the total number of such combinations . With this procedure 
the series of products that make up the £3 T s and Tl's contain no 
zero terms, and the summation is carried out in PHIBET(K) and 
TEAETA(K) over specifically defined limits. 


Subroutine 0UTHJT (IM0DE) 

This subroutine prepares the printout material. Every IPRINT 
converged solution is printed. The Fourier coefficients of the 
inplane forces, meridional transverse force, circumferential bending 
moment, twisting moment and rotations can be computed and printed 
with the solution z for the Fourier coefficients of the three 
displacements and meridional bending moment. This output material 
is converted from dimensionless form to dimensional form here. 
Provision is made to print only at stations 1, IFREQ+1, 2IFREQ+1, 
etc. This subroutine also performs the summation process for 
computing the total values of the forces, moments, displacements, 
and rotations at the NTHMAX positions around circumference prescribed 
in the input data. 


Subroutine F#LE(K) 

This subroutine prints the solution at an initial and a final 
pole . 


Subroutine THIBET (K) 

This subroutine calculates the Phis and carries out the 
multiplying and summation procedure for computing the Beta non- 
linear terms for a given meridional station K. The arrays IS, JS, 
ID, JD, IJS, MAXS, MAXD, AND MAXSY prepared in subroutine ^0DES 
are used here. 
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Subroutine TEAETA(k) 

This subroutine calculates the inplane forces and carries out 
the multiplying and summation procedure for computing the Eta non- 
linear terms for a given meridional station K. The arrays IS, JS, 
ID, JD, IJS , MAXS, MAXD, ADD MAXSY prepared in subroutine M0DES 
are used here. 


Subroutine Fj^RCE(K) 

This subroutine computes the g, (GEE), vector, equation (28), 
and the x vector, equation ( 29 a), for a given meridional station 
K. The vector GEES is the nonlinear value of g at station 1. 


Subroutine UPDATE 

This subroutine updates the storage locations of the Betas 
and Etas . It is called in subroutine XANDZ after a meridian 
station change. 


Subroutine MATIKV (A, N, B, M, DETERM, IPIV0T, INDEX, N MAX, ISCALE) 

This subroutine solves the matrix equation AX = B where A is 
a square coefficient matrix and B is a matrix of constant vectors . 

A~"^ is also obtained and the determinant of A is available. Jordan's 
method is used to reduce a matrix A to the identity matrix I through 
a succession of elementary transformations: n ]_?••• A = I. 

If these transformations are simultaneously applied to I and to a 
matrix B of constant vectors, the result is A’^ and X where AX = B. 
Each transformation is selected so that the largest element is used 
in the pivotal position. The subroutine has been compiled with a 
variable dimension statement A(NMAX, NMAX), B(NMAX, M). The 
following must be dimensioned in the calling program: IPIVOT(N MAX) , 

INDEX (N MAX, 2), A(N MAX, N MAX), B(NMAX, M) where IPIV0T and IKDEX 
are temporary storage blocks. An overflow may be caused by a 
singular matrix. The definition of the arguments of this subroutine 
are as follows : 

A = first location of a 2 -dimensional array of the A matrix. 
N = location of order of A; 

1 £ N £ N MAX 

B = first location of a 2-dimensional array of the constant 
vectors B. 
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I 


M =* location of the number of column vectors in B. 

M ~ 0 signals that the subroutine is to be used 
solely for inversion, however, in the call state- 
ment an entry corresponding to B must still be 
present. 

DETERM - Gives the value of the determinant by the 
following formula: 

DET(A) = (10 18 ) ISCALE (DETERM) 

IPIV0T - temporary storage block. 

INDEX - temporary storage block. 

N MAX - location of maximum order of A as stated in 
dimension statement of calling program. 

ISCALE - used in obtaining the value of the determinant 
by the following formula: 

DET(A) = ( 10 18 ) ISCALE (DETERM) 

At the return to the calling program A” ^ is stored at A 

and X is stored at B. 
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APPENDIX A 


CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 

The International System of Units (Si) was adopted by the Eleventh 
General Conference on Weights and Measures in i960. 

Conversion factors for the units used in this report are given in the 
following table: 


Physical quantity 

U .S . Customary Unit 

Conversion factor 
(*) 

SI Unit 
(**) 

Length 

in. 

2.54 x 10-2 

meter (m) 

Modulus of axial 
stress, elasticity 

psi 

6.895 x 10 3 

2 2 

newton/meter (N/m ) 

Temperature 

degree Fahrenheit 

K=(°F + 459.67)/l. 8 

kelvin (K) 


*Multiply value given in U.S. Customary Unit by conversion factor 
to obtain equivalent value in SI unit. 

**The prefix giga (G) is used to indicate 10^ units. 


62 



on 


APPENDIX B 


PROGRAM LISTING FOR DYNAMIC EXAMPLE 

A COMPUTER PRGGRAM FOR THE GEOMETRICALLY NONLINEAR ST AT I 
DYNAMIC ANALYSIS OF ARBITRARILY LOADED SHELLS OF REV 
CTHIS PROGRAM CONSISTS OF THE MAIN PROGRAM AND THE FOLLOWING 
C GE CM — SPECIFIES THE SHELL GEOMETRY 

C BOB— SPECIFIES THE SHELL STIFFNESSES 

C PLO AD — SP EC I FI ES THE APPLIED PRESSURE LOADS 

C TLOAD — SPECIFIES THE THERMAL LOADS 

C INITL— SPECIFIES THE INITIAL CONDITIONS 

C INLPOL 
C FNLPOL 
C MODES 
C XANDZ 
C ABC 
C PANDD 
C HJ 
C TEAETA 
C FORCE 
C UPDATE 
C MATINV 
C PMATRX 
C POLE 
C PHIBET 
C EFG 
C OUTPUT 

REAL NU, LAM, LAM2, JAY , MT , LSD18 , LSD1 N , MASS 
INTEGER SORD 

COMMON /IBL1/MNMAX/IBL2/N( 10 ) ,MNINIT/IBL3/M0,M1,M2, M3/ 
l/IBL5/IBCINL,IBCFNL/IBL6/KLL/IBL7/MNMAXO, MAXD ( 10 ) ,MAXS 
20),IS(10,10),JS(10,10),ID(10,10),JD(10,10),IJS(10)/IBL 
3IBL9/MAXM/I BLl C/ I FREQ, NTHM AX/ I BL1 1 / ICORFL , I PASS/ 1 BL 12/ 
4NC0NV/BL1/A(4,4) , BEE <4, 4) , C < 4 , 4 ) /BL3/PR ( 10) ,PX< 10) ,PT( 
5 /BL4/ P( 4 , 4, 200 ) , X ( 4 f 200 ) , ZFIM ( 4, 4 , 10 ) , Z F2 M< 4, 4, 10 ) , ZF 3 
6 ZF4M (4,4,10) 

7/BL5/TT( 10) , MT< 10 ) ,DT( 1C ) , DMT ( 10) 

9/BL7/D1 * SI 

8/BL6/Z( 4? 220) , SOE , OSE , A LOAD 
0/ BL8/R ( 200 ) , GAM ( 200 ) , OMT ( 200 ) 

1/BL9/FF$(4, 10) t ELIS (4) , GEES (4. 10) 

2/BLlO/PHIX( 1C ) , PHI T (10) ,PHI< 10) 

3/BL1I/0MXK 2C0) , PHE E, TO , T2 
COMMON/BL12/TDLI iTD 

1EL/BL13/0MEGK 4 , 4) , CAPL 1(4*4) ,0MEGL<4,4) , C APLL ( 4 T 4) , UN 
2/LAM2,LSD18 T LSDlN/BL15/NU,Ul( 10) « VI ( 10 ) « W1 ( 10 ) T V2 ( 10 ) f 
3) ,U3( 10) ,V3< 10 ) ?W3( 10) /BL16/EPS/BL17/DEL/BL18/EL1 (4 ) , E 
4TH < 6 ) 

4/ BL20/D EOMX ( 200 ) 

7/BL23/JAY(4*4) t H(4 r 4) /BL24/DL 

5(4,4,10) ,DG(4 f 4,10 J , DF ( 4, 4, 10 ) / BL 25/E ( 4, 4 ) ,F(4,4) ,G(4, 
610) ,BT3( 10) ,BXT3< 10) , BE3 ( 10 ) / BL28/ EXX3 ( 10 ) , ETT3 ( 10 ) , ET 
710) ,EX3( 10) ,ET3( 10) /BL29/BXK 10) ,BT1( 10) ,BXT1(10) ,BE1< 
8 BT2( 10 ) , BXT 2(10) , BE2 ( 10 ) / BL30/EXX 1 ( 10 ) , ETT 1 ( 10) , ETX1 ( 1 
9ETK10) , EXX2 ( 10 ) ,ETT2< 10 ) , ETX2 ( 10 ) , EXT2 ( 10 ) , EX2 ( 10 ) , ET 
0DELSQ,EXT1( 10) 

COMMON 

1 /BL34/DE E (4,4, 200 ) , DST ( 4 ,4 , 200 ) 
2/BL32/TKN,ELAST,CHAR,SIGO 
COMMON / I BL13/ ITRMAX, LSMAX 

COMMON /BL100/S0RD,TEEO/BL101/Z0(4,220) , Z2 (4 , 220 ) , Z3 ( 4 
1/BL102/DEL0AD/BL103/MASS( 200 ) 

1 /BL104/ ZDOT (4,220) 

DIMENSION S IGT ( 2 ) , SIGC ( 2 ) 

DIMENSION TITLE < 18 ) 

500 READ (5,862) TITLE 

502 READ ( 5 * 101 ) NO , SORD , I MODE , NDI MEN, NTHMAX , I FREQ , I PR INT , I 
1 »KMAX,MNMAX, M AX M, LSMAX, LCHMAX, ITRMAX, I C 
READ ( 5,102) NU,SlGO,ELAST,TKN,CHAR,TEEO 
READ( 5, 102 ) DELOAD, EPS 
R EAD ( 5 , 101 ) ITAPE 
DEL SD=DELOAD*DE LOAD 
KL=KMAX“1 
KLL=KMAX-2 
KMAX1=KMAX+1 
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I 


, ELL 
ELL 


KMAX2=KMAX+2 

AK=KL 

SIGTC 1 ) = SIGO*TKN 
S I GT ( 2 ) = SIGO/ELAST 
SIGC(l) = SIGO*CHAR/ELAST 
S IGC ( 2 ) = SIGO+TKN** 3/CHAR 
IF(NTHMAX.EQ.O) GO TO 10 
READ( 5, 105) ( TH( NTH) ,NTH=1 , NTH MAX) 

10 WRITE(6»201)N0 
WRITE(6,100J TITLE 
WRITE ( 6i 203 ) 

I F ( IBCINL.LT.O ) GO TO 11 
READ(5,103)OMEG1,CAPL1,EL1 
WRITE ( 6 , 204 ) 

WRITE (6, 205) < 0MEG1 ( 1 , J ) , J= 1 , 4) , (CAPL1! 1 , J ) ♦ J = 1 , 4 ) , ELI 
WRITE ( 6 » 206) ( 0MEG1 ( 2 , J ) , J=1 ,4) , ( CAPL1 ( 2 , J I J = 1 , 4 ! ELI 
WRITE (6, 207) < 0MEG1 (3 , J ) , J=1 , 4) , ( CAPL 1< 3 , J , J= 1 , 4) ! ELI 
WRITE! 6, 208) ( 0MEG1 (4^ J) I J=l,4) , (CAPLlUj J ) , J=ll41 ELI 
00 98 I = 1,4 
DO 98 J = 1,4 
KKLM = J/4 + 1 

0MEG1 ( I , J ) = 0MEG1 (I ,J)*SIGT(KKLM) 

98 CAPL1 ( I , J ) = CAPL 1(1, J)*SIGC(KKLM) 

14 I F C IBCFNL.LT.O) GO TO 12 

READ(5, 103)0MEGL,CAPLL, ELL 
WRITE (6,224) 

WRITE (6, 205) (0MEGH1, J) , J = 1 , 4) , ( C APLL ( 1 , J ) , J = 1 , 4 ) 
WRITE (6, 206) ( OMEGL (2, J ) , J=l, 4) , (CAPLLI 2 , J ) , J=1 , 4 ) , cul 
WRITE( 6,207) ( OMEGL ( 3 , J ) , J=1 , 4 ) , ( CAPLL ( 3 , J ) , J=1 J 4 ) , ELL 
WRITE (6, 208) ( OMEGL (4, J) ,J=1,4) , ( CAPLL ( 4, J ) ,J =1,4), ELL 
DO 99 I = 1,4 
DO 99 J = 1,4 
KKLM = J/4 + 1 

OMEGL ( I » J ) = OMEGL (I ,J)*SIGT(KKLM) 

99 CAPLL ( I , J) = CAPLL ( I , J )*S IGC ( KKLM ) 

GO TO 13 

11 WRITE (6,210) 

GO TO 14 

12 WRITE(6,211 ) 

13 WRITE ( 6 , 307 ) 

307 FORMAT ( //// ) 

IF(SORD.NE.C) GO TO 18 

it RMAX,LCHMAX,EPS 

WRITE ( 6, 213) CHAR ,TKN,ELAST, SIGO,NU 
GO TO 19 

18 WRITE(6,262)KMAX,MAXM, DELOAD, LSMAX, ITRMAX,EPS 
WRITE! 6, 263) CHAR, TKN, EL AST , S IGO , TEEO , NU 

19 IF(NTHMAX.EQ.O) GO TO 17 
WRITE (6,216) 

WRITE (6, 217 ) ( TH(NTH) ,NTH=1,NTHMAX) 

17 L AM=TKN/CHAR 
SO E = SI GO /EL A ST 
OSE=,5*SOE 
01=1, -NU 
Sl= 1 » +NU 
LAM2=LAM**2 

IF(NDIMEN.LT.l) GO TO 228 
SIG0=1. 

ELAST=1. 

TKN=1. 

CHAR=1 , 

WRITE(6,225) 

225 F0RMAT(////////////43X,34HTHE DATA IS IN NONDIMENS IONA 

14H ) 

GO TO 229 

228 WRITE (6,226) 

226F0RMAT (/ / / / / / / / / / / /45X.31HTHE DATA IS IN DIMENSIONAL F 

229 WrmoE ’ 

DO 230 M=1 , MAXM 
N(M)=0. 

PX(M)=0« 
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PT ( M)=0* 

PR ( M) «t • 

TT ( M) =0* 

MT (M) =0« 

DT < M) -0 • 

0MT<M>=0. 

MAXD(M)=0 
MAXS( M) =0 
230 MAXSYt M ) *0 
CALL GECM 

IF(SORO.NE.O) GO TO 804 
DO 86 K= 1* KMAX 
86 MASS ( K ) =0* 

WRITE ( 6 * 802 ) 

802 FORMAT ( IHlt 17Xt 15H STATION 
I6H OMEGA S 16H 


16H RADIUS 

OMEGA THETA16H 


16 

DEO 


16H RADIUS 

OMEGA THETA16H 


16 

DEO 


DO 978 K=1 * KMA X 
RKK=R (K )*CHAR 
OMXIK^OMXI OO/CHAR 
GAMK=GAM< K ) /CHAR 
GMTK=OMT <K)/CHAR 
D£QMXK=DEOMX (K ) / i CHAR* CHAR ) 

978 WR ITE { 6 t 803 ) K t RKK ,GA MK , OMXI K t OMTK* DEOMXK 

803 FORMAT < 2GX t 13 » 9X,5E16.4) 

GO TO 805 

804 WRITE i 6 » 810 ) 

810 FORMAT ( 1H1 » 5X , 1 5H STATION 

1 16H OMEGA S 16H 

2 MASS //) 

TEED-TEEO 

IF(NDIMEN#EQ#1) TEED^l. 

DO 979 K=1 » KMAX 

RKK=R<K)*CHAR 

OMX IIOOMXI <K)/CHAR 

GAMK=GAM (K ) /CHAR 

GMTK=OMT <K) /CHAR 

DEOMXK=DEOMX(K ) / < CH AR*CHAR > 

AMSS~MA$S(K)*TEED**2*ELAST*TKN/CHAR**2 

979 WRITE ( 6 » 81 3 ) K , RKK ,GA MK , OMXI K , OMTK, DEOMXK , AMS S 
813 FORMAT ( 8X,I3 f 9X t 6E16.4> 

805 M0=0 
M1=0 
M2=0 
M3=0 

ABN-CHAR/S I GO/TKN 
ZN=S IGO*TKN 
WRITE(6 t 112) 

112 FORMAT <////17X*12H STATION 20H 8 STIFFNESS 2 

1 1 FFNESS 20H B PRIME 20H D PRIME 

DO 888 K=l f KMAX 
CALL BDB(K, BtDBtD, DD) 

BST=ELAST*TKN 
Z$T=ELAST*TKN**3 
B=B*B$T 
D = D*ZST 
DB=DB/CHAR*B$T 
DD=DD/CHAR*ZST 
888 WRITE ( 6 t 71 ) K t B,D,DB f DD 
71 FORMAT ( 20X » 1 3 , 4X ? 4E 20* 6 > 

CALL PLOAD(l) 

CALL TLOAO(l) 

I F ( SORD* NE • 0 ) GO TO 891 
DO 889 M=1 T MNMAX 
WRITE 16,113) N(M) 

113 F0RMAT(///25Xt44HPRESSURE AND TEMPERATURE COEFFICIENTS 
1 FOLLOW//) 

WRITE (6*114) 

114 FORMAT <5Xf7HSTATIQNt3Xtl5H PR 1 5H PX 

1 PT 15H TT 15H MT 15H 

25H DMT //) 

DO 890 K=1 t KMAX 
CALL PLOAD(K) 
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CALL TLOAD(K) 

PRM=PR ( M ) / ABN 
PTM=PT(M)/ABN 
PXM=PX<M)/ABN 
TTM=TT ( M)*ZN 

ENTM=MT(M)/CHAR*ZN*TKN*TKN 
DTM=DT ( M )/CHAR*ZN 

DMTM=DMT { M ) *ZN*TKN*TKN/ ( CHAR*CHAR ) 

890 WRITE (6,115) K , PRM, PXM , PTM , TTM , EMTM , DTM , D MTM 
115 FORMAT ( 6X, I 3, 7X , 7E 1 5* 4 ) 

889 CONTINUE 

891 CONTINUE 
DELSQ=DEL**2 
TDLI=. 5/DEL 
TDEL=2.*DEL 
MNINIT=1 
MNMAXO=MNMAX 
DO 20 1=1,4 
DO 20 J = 1 , 4 
UNIT ( I , J )=0« 

IF(I.EQ.J) UNIT( I , J ) =1 • 

20 CONTINUE 

NMAX=MAXM*KMAX2 
DO 22 K= 1 , NMAX 
DO 22 1=1,4 
ZDOT ( I * K ) =0« 

Z2( I,K)=0. 

Z3( I,K)=0. 

Z0(I,K)=0. 

22 Z( I » K ) =0 • 

ALOAD=DELOAD 

IF(IC.EQ.O) GO TO 834 

I F( IT APE«GT • 1 ) WRI TE( 6, 280) 

280 FORMAT ( ////5X,28HTHIS IS A RESTARTED SOLUTION//) 
IF(ITAPE.GT.l) GO TO 834 
CALL INITL 
ACO=CHAR*SIGO/ELAST 
ACM=S I GO*TKN** 3/CHAR 
DO 830 M=1 , MNM AX 
MM=(M-1 )*KMAX2 
WR I TE ( 6 , 126 ) N ( M ) 

126 FORMAT (////5X,29HTHE INITIAL CONDITIONS FOR N=I3,8H F 
WRI TE ( 6 , 127 ) 

127 FORMAT ( 19X,7HSTATI ON, 3X,20H U 20H 

1 20H W 20 H M S //) 

DO 831 K=2 , KMA XI 
MK=K+MM 

TU=ACO*ZO( 1 , MK ) 

TV=AC0*Z0(2,MK) 

TW=AC0*Z0(3,MK ) 

TM=ACM*Z0(4,MK) 

WRITE(6,71) KK,TU,TV,TW,TM 

831 CONTINUE 
WRITE! 6» 129) 

129 F0RMAT(///19X»7HST ATI ON * 3X » 20H U DOT 2CH 

1 2 OH W DOT 20H M S DOT / 

DO 832 K=2 * KMA X 1 
ACD=CHAR*SIGO/( ELAST*TEEO) 

AMD=SIG0#TKN**3/ (CHAR*TEEO) 

MK=K+MM 

TU=ACD* ZDOT ( 1 , MK ) 

TV=ACD*ZD0T(2,MK) 

TW=ACD* ZDOT < 3, MK ) 

TM= AMD* ZDOT (4» MK ) 

WRITE ( 6 , 71 ) KK, TU, TV, TW,TM 
DO 833 1=1,4 

Z( I,MK)=ZO( I , MK ) +ZDOT ( I , MK ) *DELOAD 
Z2( I , MK ) =Z0 ( I , MK ) — Z DOT ( I,MK) *DELOAD 
833 Z3II ,MK)=Z0II,MK)-2,*ZD0T< I,MK)*DELOAD 

832 CONTINUE 
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31 Z(It IK)-ZO( IvlK) 

GO TO 62 

361 WRITE ( 6 » 271 ) 

WRITE RESTART DATA ON TAPE UNIT 8 
IF( ITAPE.NE.C) REWIND 8 

IFUTAPE.EQ.l.OR. I TAPE. EQ. 3) WRITE (8) ! t Z ( I , JP ) , ZO { I , 

1 ,Z3( I,JP) ,1=1,4) , JP=1,NMAX) 

GO TO 500 

365 WRITE ( 6,266) ITRMA X , LSTEP , NO 
GO TO 500 
862 FORMAT (18A4) 

100 FORMAT (38X,18A4,///) 

101 FORMAT (1615) 

102 FORMAT ( 6E12.3 ) 

103 FORMAT !4E16.8) 

105 FORMAT ( 6E12. 3) 

201 FORMAT ( 1H1 » 48X , 16H — PROBLEM NUMBER, 14, 2H — ,///) 

202 FORMAT ( 1 X12 A6/ / / ) 

203 FORMAT !49X,21H — INPUT DATA RECORD — ,//34H THE BOUN 

IONS ARE: // ) 

204 FORMAT ! 29H AT THE INITIAL EDGE//2X, 40H- 


-LAMBDA BAR- 


N S 


< ,4E10.3,12H) U 
( ,4E10.3,12H) V 


< »4El0.3f 12H) 
( ,4E10.3,12H) 


1A BAR- 1 5Xt 40H- 

2 — 12X 1 1 OH EL /) 

205 FORMAT ! 1 X , 1H ( ,4E10.3,15H) 

1 E10.3) 

206 FORMAT (1X»1H( ,4E10.3,15H) N ST 

1 E10.3) 

207 FORMAT ( IX, 1H! ,4E10.3, 15H) Q S 

1 E10.3) 

208 FORMAT! IX, 1H! ,4E10.3,15H) PHI S 

1 E10 • 3 ) 

210 FORMAT ( 39H THE SHELL HAS AN INITIAL POLE/) 

211 FORMAT ( 1H0 , 36H THE SHELL HAS A FINAL POLE//) 

212 FORMAT ( 5X , 40HNUMBER OF STATIONS 1 

1ER OF MODES 13/ 5 X , 40H INC REMENT 

20R F6.3/5X,40HMAXIMUM NUMBER OF LOAD S 

3 1 3 / 5X , 40 H MAX I MUM NUMBER OF ITERATIONS 1 

4MUM NUMBER OF LOAD FACTOR CHANGES 1 3/ 5X , 40HC0NV ERGEN 

5 F6.3//) 

213 FORMAT! 5X ,33HC HARACTE RI ST I C SHELL DIMENSION E12.4/ 5 

ICE THICKNESS E 1 2.4/5X , 33HREFERENCE ELASTI 

2 E12.4/5X, 3 3 HR E FERENC E STRESS -E12. 

3 SON 1 S RATIO E12.4//) 

216 FORMAT! 5X, 62HC I RCUMF ERENTI AL COORDINATES FOR PRINT RE 
1 IANS, ARE:/) 

217 FORMAT (10X,E16.6,5(1H» ,E16.6,1H )) 

220 FORMAT! 1H1,80H THE MAXIMUM NUMBER OF LOAD CHA 

1EN MADE. END PROBLEM NUMBER 14 ) 

221 FORMAT! 1H1,79H THE MAXIMUM NUMBER OF LOAD STE 

1 TAKEN. END PROBLEM NUMBER 1 4 ) 

222 FORMAT ! 1 HI , 1 19H THE SOLUTION DID NOT CONVERGE 

1AXIMUM NUMBER OF ITERATIONS. THE LOAD FACTOR HAS BEEN 
25.) 

223 FORMAT (1H1,69H THE SOLUTION DID NOT CONVERGE 

IT LOAD INCREMENT. /11X,71HL00K FOR AN ERROR IN THE INPU 
2RY A SMALLER VALUE FOR DELOAD.) 

224 FORMAT! 1H0,26H AT THE FINAL EDGE//2X,40H 

IGA BAR 1 5X, 40H LAMBDA BAP 

2 >12X,10H EL /) 

262 FORMAT! 5X.40HNUMBER OF STATIONS 1 

1ER OF MODES 1 3/ 5X , 40H INCREMENT 

2— E9. 3 / 5X , 40HM AX I MUM NUMBER OF TIME S 

3 I5/5X, 40HMAXI MUM NUMBER OF ITERATIONS 1 

4ERGENCE CRITERION F6.3//) 

263 FORMAT (5X,33HCHARACTERISTIC SHELL DIMENSION E12.4/ 5 

ICE THICKNESS E12.4/5X, 33HREFERENCE ELASTI 

2 E12.4/5X, 33HRE FERENC E STRESS E12. 

3RENCE TIME El 2. 4/ 5X , 33HP0I S SON 1 S RAT 

4 El 2 .4/ / ) 

266 FORMAT! 1H1,35H THE SOLUTION DID NOT CONVERGE INI3,24 
1 AT TIME STEP I 5 , 21H. END PROBLEM NUMBE R I 4 , 1 H. ) 

271 FORMAT! 1H1,79H THE MAXIMUM NUMBER OF TIME STE 
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1 TAKEN 
END 


END PROBLEM NUMBER 14 ) 


SUBROUTINE GEOM 
REAL MASS 
COMMON 

1 / 1 BL4/KMAX » KL 

0 /BL 8/R ( 200 ) , GAM (200) ,0MT<200) 
3/BL11/OMXH200) ,PHEE,T0,T2 
4/ BL17/DE L 
4/BL20/DE0MX ( 200 ) 
6/BL32/TKN,ELAST,CHAR,SIG0 
I/BL102/DELOAD/BL103/MASS(200) 

C GEOMETRY DATA 
AKX=KMAX-1 
DE L = 1 • / AKX 

THET=ARSIN( 2.2801/CHAR) 

DO 11 K=1 , KMAX 
AK=K 

R ( K ) = ( 7 • 9499+ ( AK-1 . ) * ( 2. 280 1 ) / AKX ) /CHAR 
GAM( K)=( 2.2801 /CHAR)/R(K) 

OMX I ( K ) =0. 

DEOMX ( K ) =0 . 

OMT(K)=COS(THET)/R(K) 

MASS ( K ) = 1. 

11 CONTINUE 
RETURN 
END 


SUBROUTINE BDB ( K , B , DB , D , DD ) 

REAL NU 
COMMON 

1/BL32/TKN,E LAST, CHAR, SIG0/BL15 

2/NU,Ul( 10) , VI ( 10) ,W1< 10) ,V2< 10) ,U2( 10) ,W2(10) ,U3< 10 ) ,V 
3 / BL17/DEL 
C STIFFNESS DATA 
B=l. 089082 
D=. 09075683 
DB=0. 

DD=0. 

RETURN 

END 


SUBROUTINE PLOAD(K) 

COMMON 

1/BL32/TKN,ELAST,CHAR,SIG0/IBL2/NN( 10) ,MN I NIT 

8 / BL6/ Z (4,220) , S OE , 0 SE , ALO AD 

4/IBL1/MNMAX 

5/BL3/PR( 10) ,PX( 10) t PTC 10) 

COMMON/ IBL4/KMAX,KL 

C0MMON/BL8/R(200) ,GAM( 200) ,OMT( 200) 

1 / I BL8/L STEP , ITR 
1/BL102/DEL0AD/BL 103 /MASS (200) 

RETURN 

END 


SUBROUTINE TLOAD(K) 

REAL NU 
COMMON 

1/IBL1/MNMAX/IBL2/NN(10) , MNINI T 
2/BL5/TT ( 10) ,EMT( 10) ,DT( 10) , DMT ( 10) 

3/BL32/TKN, ELAST , CH AR , S I GO 
4/BL6/Z( 4,220) , SOE , 0 SE , ALOAD/ B LI 5/ 

5NU.UK 10) , VI ( 10 ) ,W1 ( 10 ) ,V2(1C ) ,U2( 10) ,W2< 10) , U3 ( 10) ,V3 
1/IBL8/LSTEP.ITR 
RETURN 
END 
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SUBROUTINE INI TL 

COMMON /BL101/ ZO ( 4, 220 ) , Z2 < 4, 220 ) , Z3< 4, 220) , DEL SD 

1 /BL104/ZD0T(4t22CI) /BL6/Z (4» 220 ) »SOE»OSE» ALOAD 
1/IBL2/NNI 10 ) iMNI NIT 

2 /I BL1/MNMAX/I BL9/MAXM/IBL12/KMAXltKMAX2 » NCONV/I BL4 /KM 
3/BL32/TKNtELAST . CHAR* S I G0/BL1 OO/SORD »TEEO 

INITIAL CONDITIONS DATA 
NN ( 1 ) =0 
NN ( 2 ) =1 
NN( 3)=2 
NN( 4 ) =4 
PI=3. 14159 
DO 2 M=1.MAXM 
IF(M.EQ.l) VEL=-444. 08/PI 
IF(M.EQ.2) VEL=-444.08/2. 

IFIM.EQ.3) VEL=-444.08*2./(3.*PI) 

IFIM.EQ.4) V EL =444. 08*2. /(15.*PI ) 

DO 2 K=2,KL 
I=K+ l+( M— 1 ) *KMAX2 

2 ZDOT ( 3» I ) = VEL*ELAST*TEEO/ ( CHAR* S I GO) * 10 • 

RETURN 

END 


SUBROUTINE INLPOL 
COMMON 

1 / I BL1/MNMAX 
l/IBL3/M0,Ml,M2tM3 
2/IBL12/KMAXl»KMAX2» NCONV 
8 /BL6/Z ( 4 » 220 ) , SOE , OSE , ALOAD 
4/BL7/D1 » SI 

3/BL11/0MXK 200) tPHEE,T0,T2 
6/BL 17/DEL 

7/BL29/BXK 10) , BT1 ( 10) , BXT 1 < 10 ) , BE 1 < 10 ) , BX2 ( 1C ) , BT2( 10) 
8BE2I 10) 

9/BL30/EXX1I 10) .ETTK10) ,ETX1(10) , EX1 ( 10 ) , ET1 ( 10 ) , EXX2 ( 
0,ETX2(10) .EXT2( 10) , EX2 ( 10 ) , ET2 ( 10 ) /BL 31 /DELSQ t EXT1 ( 10 ) 
1/BL5/TT (10) ,EMT(10) , DT ( 10 ) , DMT ( 10 ) 

2/IBL13/ITRMAX,L$MAX 
DO 1 MN=1 » MNMAX 
BX1 ( MN ) =0. 

BT1 ( MN ) =0 . 

BXT1< MN )=0. 

BE1 ( MN ) =0. 

EX1 ( MN )=0. 

ET1 ( MN ) =0. 

ETX1 ( MN ) =0. 

L EXX1 ( MN ) =0. 

IF(Ml.EQ.O) RETURN 

I2=2+(M1-1)*KMAX2 

13=12+1 

14=13+1 

PHEE=(1.5*Z(3, 12 )-2.*Z<3, I3)+.5*Z(3, 14) ) /DEL+OMXI ( 1 ) *Z 

BET=.5*PHEE**2 

I F ( ITRMAX. EQ.l ) BET =0 . 

T2=0. 

IF ( M2. EQ.O ) GO TO 2 
CALL BDB(1»B»DB»D» DD) 

I2=2+(M2-1)*KMAX2 

13=12+1 

14=13+1 

T2=B*D1*( (-1.5*Z(1,I2)+2.*Z(1,I3)-.5*Z(1, 14) )/DEL+.5*S 

Q1=.5*PHEE*T2 

BX1 ( M2 ) =BET 

BT1 ( M2 ) =-BET 

BXT1IM2 )=-BET 

ETX1 ( Ml ) =Q1 

I F ( M3. EQ.O ) GO TO 2 
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EXX1 ( M3 ) =Q1 
ETX1 ( M3 ) =— Q1 

2 T0=0. 

IF(MO-EQ.O) GO TO 3 
BX1 ( MO) =BET 
BT1 ( MO ) =BET 
CALL BDB(l t B*DB,D,DD) 

CALL TLOADU) 

I2=2+( MO— 1 ) #KM AX2 
13=12+1 

TO=B*Sl*( (-1.5*Z(lt I2)+2»*Z<1» I3)-.5*Z< 1 , I A) J /OEL+OMXI 
1 + . 5*SOE*BET ) -TT ( MO > *ALOAD 

3 EXX1(M1)=PHEE*(T0+.5*T2) 

RETURN 

END 


SUBROUTINE FNLPOL 
COMMON 
1/IBL1/MNMAX 
2/IBL3/M0,Ml,M2,M3 
3/IBL4/KMAX,KL 
4/IBL12/KMAXl,KMAX2,NC0NV 
8/BL6/Z(4i220) t SOE , OSE , ALOAD 
6/BL7/D1 ,S1 

3/BLU/OMXK 200) ,PHEE,T0,T2 
8/BL17/DEL 

9/BL27/BX3( 10) »BT3( 1 0 ) ,BXT3( 10) , BE 3 ( 10) 

0/BL28/EXX3< 10 i , ETT3 < 10 ) , ETX3( 10),EXT3(10) , EX3 ( 10 ) , E T3( 
1/BL5/TT( 10) ,EMT( 10) »DT< 10) , DMT ( 10) 

2/IBL13/ITRMAX, LSMAX 
DO 1 MN = l,Mt\MAX 
BX3 ( MN ) =0. 

BT3 ( MN ) =0. 

BXT3( MN ) =0. 

BE3 ( MN ) =0. 

EX3 (MN)=0. 

ET3 ( MN ) =0. 

ETX3(MN)=0. 

L EXX3(MN)=0. 

CALL BDB(KMAX,B,DB,D,DB) 

IF(Ml.EQ.O) RETURN 
KM=KMAX1+ ( Ml-1 ) *KMAX2 
KM1=KM— 1 
KM2~K M 2 

PHEE=-< 1.5*Z(3 , KM)-2.*Z< 3*KM1 )+. 5*Z(3 ,KM2 ) )/DEL+OMXI (K 

BET=.5*PHEE**2 

IF( ITRMAX.EO.l ) BET=0. 

T2=0. 

I F ( M2 « EQ.O ) GO TO 2 
KM=KMAX1+(M2-1 )*KMAX2 
KM1=KM-1 

2 

T2=B*D1*< (1.5*Z(1,KM1-2.*Z<1,KM1)+.5*Z< 1,KM2) )/DEL+.5* 
Q1=.5*PHEE*T2 
BX3 ( M2 ) =BET 
BT3 ( M2 ) =— BET 
BXT3 ( M2 ) = BET 
ETX3( Ml )=Q1 
I F ( M3 . EQ.O ) GO TO 2 
EXX3(M3 ) =Q1 
ETX3( M3 ) =— Q1 
2 T0=0. 

IFCMO.EQ.O) GO TO 3 
CALL TLOAD(KMAX) 

KM=KMAX1+(M0-1 )*KMAX2 

KM1=KM- 1 

KM2=KM-2 

BX3 ( MO ) =BET 

RT3 f MO ) —BET 

T0=B*S1*( <1.5*Z(1,KM)-2.*Z< 1,KM1)+.5*Z< 1,KM2) J/DEL+OMX 
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1Z(3,KM)+.5#S0E*BET) -TT ( MO ) *ALOAD 
3 EXX3C Ml ) =PHEE* ( TO+« 5*T2 ) 

RETURN 

END 


SUBROUTINE MODES 
COMMON 

1/IBL1/MNMAX 

2/ I BL2/N ( 10 ) tMNINIT 

3/1 8L7/MNMAX0 , MAXD ( 10) * MAXS ( 10) ,MAXSY<10) , IS (10,10) , JSC 
4«10)«JD(10 f 10)f IJS(IO) 

5/IBL9/MAXM 
6/IBL11/ I CORF L, IPASS 
IF(MAXM.EQ.l) RETURN 
IF(MNINIT.GT.MAXM) RETURN 
DO 1 MN=1, MNMAXO 
NMN=N( MN ) 

NNS=MN 

IF(MNINIT.GT.MN) NNS=MN INI T 
DO 1 MM=NNS, MNMAXO 
NMM=N( MM ) 

NTEST= I ABS ( NMN-NMM ) 

DO 2 MMFT=1 ,MNMAX 
IFCNTEST.EQ.N(MMFT) ) GO TO 10 
2 CONTINUE 

I F ( ICORFL.EQ.l ) GO TO 1 
MNMAX=MNMAX+1 
N( MNMAX ) =NTEST 
MMFT=MNMAX 

IF(MNMAX.EQ.MAXM) IC0RFL=1 

10 I F ( NMN-NMM ) 11,1,12 

11 LOCD=MAXD( MMFT ) +1 
MAXD(MMFT)=LOCD 

I D ( LOCD, MMFT )= MM 
JD( LOCO , MMFT) =MN 
GO TO 1 

12 L0CD=MAX0(MMFT)+1 
MAXD( MMFT) =LOCD 
ID( LOCD , MMFT )=MN 
JD( LOCD , MMFT )=MM 

1 CONTINUE 

DO 301 MN=1, MNMAXO 
NMN=N( MN) 

NNS=MN 

IF(MNINIT.GT.MN) NNS=MNI NIT 
DO 301 MM=NNS, MNMAXO 
NMM=N( MM ) 

NTEST=NMN+NMM 
DO 302 MMFT=1, MNMAX 
I F ( NTEST* EQ* N( MMFT ) ) GO TO 310 
302 CONTINUE 

IF( ICORFL.EQ.l) GO TO 301 
IF( MNMAX.GE.MAXM) GO TO 301 
MNMAX=MNMAX+1 
N (MNMAX ) =NTE ST 
MMFT=MNMAX 

IF ( MNMAX.GE.MAXM) IC0RFL=1 
310 IF(NMN.EQ.NMM) GO TO 360 
LOCS=MAXS ( MMFT ) +1 
MAXS ( MMFT) =LOC S 
ISC LOCS , MMFT ) = MN 
JS ( LOCS , MMFT ) =MM 
GO TO 301 

360 IF(NMN.EQ.O) GO TO 301 
MAXSY ( MMFT ) =1 
I JS ( MMFT)=MN 
301 CONTINUE 

MNINIT=MNMAX0+1 

IFC ICORFL.GT.O) I PASS=I PASS+1 

IF( IPASS. LT. 2. AND. MNINIT.LE. MNMAX) CALL PMATRX 
RETURN 
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END 


SUBROUTINE XANDZ 
REAL NU , LAM2tMT * MASS 

COMMON/ BL5/TT( 10) ,MT(10),DT< 10 > » DMT { 10) 

COMMON 
1/ I BL1/MNMAX 
2/1 BL3/M0 1 Ml * M2 ? M3 
3/1 BL4/KMAX ? KL 
4/IBL5/IBCINL, IBCFNL 
5/IBL6/KLL 

3/IBL7/MNMAX0,MAXD( 10) . MAXS ( 10) , MAXSY ( 10 ) , IS ( 10,10 ), JS ( 
4,10) ,JD(10,10) ,IJS( 10) 

COMMON/IBL8/LSTEP t itr 
8/ IBL12/KMAX1 » KMAX2 » NCONV 
1/IBL13/ITRMAX, LSMAX 
9/BLl/A<4»4) ,BEE(4,4) ,C(4,4) 

2/BL3/PR( 10 ) , PX ( 10) f PT ( 10) 

5/BL4/P(4»4»200 ) , X ( 4 , 200 ) , ZF 1M ( 4 , 4, 1C ) , ZF2 M( 4, 4, 10 ) , ZF3 
1 ZF4M (4,4,10) 

8/BL6/ZI 4,220 ) ,SOE,OSE, ALOAD 
3/BL7/D1 ,S1 

0 /BL8/R ( 200 ) ,GAM( 200) ,OMT< 200) 

5/BL9/FFS(4,10) , ELIS (4) t GEES (4, 10) 
6/BL14/LAM2,LSD18,LSDlN 

7/ BL1 5/NU , U1 ( 10) , Vl< 10) ,W1< 10) ,V2( 10) ,U2< 10) ,W2(1C) ,U3< 
8 W3 ( 10 ) 

COMMON 

1/BL16/EPS 

2/BL27/BX3U0) ,BT3< 10) ,BXT3(10) ,BE3< 10) 

3/BL28/EXX3 ( 10) ,ETT3(10) ,ETX3( 10) ,EXT3< 10) ,EX3(10) ,ET3( 
4/BL29/BX1C 1C) ,BT1( 10) ,BXT1(10) ,BE1( 10) ,BX2( 10) ,BT2( 10) 
5BE 2 ( 10 ) 

6/BL30/E XXI ( 10) * ETTl(lO) « ETX1 ( 10 ) , EX1 < 10 ) , ET1 ( 10) , EXX2 < 
7 y ETX2 ( 1 0 ) • E XT2 ( 10) , EX2( 10) ,ET2< 10) 

8/BL31/DELSQ,EXTl ( 10 > 

9/BL18/ELK4), ELL (4) 

COMMON / BL1 00/ SORD ,TEEO/BL101/Z0(4,220) , Z 2 ( 4 , 220 ) , Z 3 < 4 
1 /BL104/ZD0T<4,220) /B L102/DE LOAD 
1/BL103/MASS(200) 

DIMENSION ELLS (4) ,FLS (4) ,ZT(4) , IPIV0T(4) , INDEX (4, 2) 

1, CL0(4,4) tCLl( 4,4) , CL 2 ( 4 » 4 ) 

2, TZMAX(4,10),ZDD(4) 

EQUIVALENCE ( C LC ( 1 ) , ZF 1 M (1) ) , < CL1 < 1 ) , ZF 2 M ( 1 ) ) , ( CL2 ( 1 ) , 

DO 201 1=1,4 

DO 201 M=1 , MNM AX 

MJ=1+(M— 1)#KMAX2 

TZMAX( I f M) = ABS (Z< I , M J ) ) 

DO 201 K=2 , KMAX2 
KM=K+(M-1)*KMAX2 
AZTST=ABS(Z( I,KM) ) 

IF(AZTST.GT.TZMAX( I,M) ) T ZM AX ( I , M ) = AZT ST 
201 CONTINUE 
NC0NV=1 

IF ( ITRMAX. EQ. 1 ) GO TO 66 
DO 1 M= 1 , MNMAXO 
I=1+(KMAX+2)*(M-1) 

U1(M)=Z( 1,1 ) 

VI ( M) =Z ( 2 1 I ) 

W1(M)=Z<3,I ) 

11=1+1 

U2(M)=Z< 1,11) 

V2 ( M ) =Z (2,11) 

1 W2(M)=Z(3, ID 
IFUBCINL.LT. 0) GO TO 100 
CALL PH I BET ( 1 ) 

DO 2 M= 1 , MNMAX 
BX1(M)=BX3(M) 

BT1 ( M ) =BT3( M ) 

BXT1 (M)=BXT3(M ) 

2 BE1 ( M ) =BE3 ( M ) 
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3 

102 

4 


5 

66 


67 

8 

9 

20 


10 


11 

12 
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) 

GO TO 

20 

DB 

,D, DD) 


) 

GO TO 

67 

>*ALOAD+OSE* 

B1 

* D1 



CALL TEAETA ( 1 ) 

DO 3 M=1 1 MNMAX 
EXX1(M)=EXX3(M) 

ETT1 ( M ) = ETT3( M ) 

ETX1 ( M) =ETX3 ( M ) 

EXT1 ( M) =EXT3 ( M ) 

EX1 ( M) = EX3 ( M ) 

ET1 ( M ) = ET3 ( M ) 

CALL PHI BET ( 2 ) 

DO 4 M=l, MNMAX 
BX2 ( M ) = BX3 ( M ) 

BT2(M)=BT3(M) 

BXT2 ( M ) =BXT3 ( M ) 

BE2 ( M) = BE3 ( M ) 

CALL TEAETA( 2 ) 

DO 5 M=l, MNMAX 
EXX2 (M)=EXX3(M> 

ETT2( M)=ETT3 ( M ) 

ETX2( M) =ETX3( M ) 

EXT2(M)=EXT3(M) 

EX2 ( M ) = EX3 ( M ) 

ET2<M)=ET3<M) 

CALL PHIBET ( 3) 

CALL TEAETA(3 ) 

CONTINUE 
IF( IBCINL«LT#0) 

CALL BDB < 1 , B1 , 

GAM1=GAM( 1 ) 

CALL TLOAD(l) 

DO 8 M=l, MNMAX 

FFS( 1?M )=-TT(mUaL0AD+0SE*(BX1(M)+BE1(M)+NU*(BT1( M)+BE 
FFS(2,M)=0SE*( B1 *D1 * BXT1 ( M ) + EX1 < M ) +ET1 ( M ) ) 

FFS(3,M)=LAM2*GAM1*D1*MT(M)*AL0AD-(EXX1< MJ+ETXl (M) )*S0 
GO TO 8 

FFS ( 1 » M ) =-TT ( M) * ALOAD 
FFS(2,M)=0. 

FFS ( 3 , M ) =LAM2*G AM1*D1*MT < M ) *ALOAD 
FFS<4,M)=0. 

DO 9 1=1.4 

EL1S<I) =AL0AD*EL1 ( I ) 

CALL FORCE ( 1 ) 

CALL FORCE ( 2 ) 

DO 10 K=3» KLL 
KP=K+1 

IF( ITRMAX.EQ.l ) 

CALL UPDATE 
CALL PHIBET ( KP ) 

CALL TEAETA(KP) 

CALL FORCE ( K) 

I F ( ITRMAX.NE.l ) 

IF(IBCFNL.LT.O) 

IF ( ITRMAX# EQ* 1 ) 

CALL PHIBET ( KMAX ) 

CALL TEAETA(KMAX) 

CALL FORCE ( KL ) 

CALL FORCE(KMAX) 

DO 12 1=1,4 

ELLS ( I ) =ALOAD*ELL< I ) 

CALL BDB (KMAX, BL,DB,D, DD) 

GAM L=GAM( KMAX) 

FLS ( 4) =0* 

CALL TLOAD(KMAX) 

DO 14 M=1 T MNMAX 
I F ( ITRMAX# EQ« 1 ) GO TO 68 

FLS< 1 )=-TT(M)*AL0AD+0SE*(BX3(M)+BE3<M)+NU*( BT3{M)+BE3( 
FLS ( 2 ) =0SE*(BL*D1*BXT3(M)+EX3(M)+ET3(M) ) 

FLS(3) = LAM2*GAML*Dl*MT(M)*ALOAD-< EXX3 ( M ) +ETX3 ( M ) )*SOE 
GO TO 69 

FLS ( 1 )=— TT ( M )*ALOAD 
FLS ( 2 ) =0 • 

FLS(3)=LAM2*GAML*D1*MT(M)*AL0AD 


GO TO 10 


CALL UPDATE 
GO TO 120 
" TO 11 


GO 
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69 CONTINUE 

IK=KL+KMAX*(M-1 ) 

I J=KMAX*M 
L=M*KMAX2 
DO 14 1=1,4 
SUMZ=0. 

DO 15 J = 1 , 4 

15 SUMZ=$UMZ+ZF1M< I, J,M)*ELLS< J)+ZF2M< I , J , M ) *X( J , I J ) +ZF3M 
IX ( J, IK)+ZF4M( I , J,M)*FLS< J) 

14 Z ( I , L ) = $UMZ 
LS=1 

150 DO 16 M=l, MNMAX 
DO 16 L=LS , KMAX 
K=KMAX2-L 
KPX=K— 1 
KZ=K+1 

IJ=KPX+(M-1)*KMAX 
JK=KZ+(M-1)*KMAX2 
KK= JK— 1 
DO IT 1*1,4 
SUMZ=0. 

DO 18 J = 1 , 4 

18 SUMZ=SUMZ-P( I, J, IJ)*Z( J, JK) 

SUMZ=SUMZ+X ( 1 1 I J ) 

ASUMZ=ABS( SUMZ ) 

I F ( ASUMZ.GT.l.E+15) ITR= I TRMAX 

I F ( NCON V • NE • 1 .OR. ASUMZ .LT. 1.E-D5) GO TO 17 
DEL Z= A B S ( Z ( I,KK)-SUMZ) 

ZTEST=EPS#TZMAX( I , M) 

IF(DELZ.GT.ZTEST) NC0NV=0 
17 Z ( I , KK ) = SUMZ 

16 CONTINUE 

I F ( IBCINL.LT. 0) GO TO 30 
DO 25 M=l, MNMAX 
CALL EFG ( 1 , M ) 

CALL ABC 

IJ=2+(M-1)*KMAX2 
I Jl=I J+l 
I J2= I J- 1 
DO 21 1=1,4 
SUMZ=0. 

DO 22 J= 1 , 4 

22 SUMZ= SUM Z- A ( I, J )*Z( J, I J 1 ) -BE E ( I , J ) *Z ( J , I J ) 

21 ZT ( I ) =SU MZ+GEE S ( I , M ) 

CALL MAT I NV ( C , 4 , ZT , 1 , DETE RM , I P I VOT , I NDE X , 4 , I SC ALE ) 

DO 23 1=1,4 

23 Z ( I , I J2 ) =ZT ( I ) 

25 CONTINUE 

30 RETURN 

100 CALL I NL POL 

DO 101 M=1 « MNM A XO 
U1 ( M)=U2 (Mi 
Vl< M ) =V2 ( M ) 

W 1 ( M ) =W2 ( M ) 

I J=3+KMAX2*(M-1 ) 

U2 ( M ) =Z ( 1 , I J ) 

V2 ( M ) =Z ( 2 , IJ) 

101 W2 ( M ) =Z ( 3 , I J ) 

GO TO 102 

120 I F < ITRMAX.NE.l ) CALL FNLPOL 
CALL FORCE! KL) 

I F ( M2 * EQ.O ) GO TO 122 
L=KL+(M2-1)*KMAX 
L 1=KMAX 1 +( M2-1 )*KMAX2 
DO 130 1=1,4 
SUM=0. 

DO 131 J=1 , 4 

131 SUM=SUM+CL2(I, J)*X< J,L) 

ASUMZ= ABS (SUM) 

I F ( NCON V. NE • 1 .OR. ASUMZ.LT. l.E-05) GO TO 130 
DELZ=ABS ( Z ( I, LI )-SUM) 

ZTEST=E PS*T ZMAX ( I, M2) 
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IF(DELZ.GT.ZTEST) NC0NV=0 
130 Z ( I »L1 ) = SUM 

122 IF(Ml.EQ.O) GO TO 123 
L=KL+( M1-1)*KMAX 
Ll=KMAXl+( Ml— 1 )*KMAX2 
DO 132 1=1,4 

SUM=0. 

DO 133 J=1 , 4 

133 SUM=SUM+CL1(I, J)*X( J,L) 

ASUMZ=ABS(SUM) 

IF(NCONV.NE.l .OR. ASUMZ .LT. l.E-05) GO TO 132 
DELZ=ABS ( Z ( I » L 1 )-$UM) 

ZTE$T=EPS*TZMAX( I , Ml ) 

IF(DELZ.GT.ZTEST) NC0NV=0 
132 Z( 1 1 LI ) = SUM 

123 IF(MO.EQ.O) GO TO 124 
L=KL+(M0-1)*KMAX 
L1=KMAX 1 +( MO— 1 )*KMAX2 
DO 134 1=1,4. 

SUM=0. 

DO 135 J = 1 , 4 

135 SUM=SUM+CLO ( I , J ) *X f J , L ) 

ASUMZ=ABS( SUM) 

IF ( NCONV.NE. 1 .OR. ASUMZ. LT . 1 . E-Q 6 ) GO TO 134 
DELZ=ABS ( Z ( I , LI ) -SUM ) 

ZTE$T=EPS*TZMAX( I , MO) 

IF(DELZ.GT.ZTEST) NCONV=0 

134 Z ( I , LI ) =SUM 

124 LS=2 

GO TO 150 
END 


SUBROUTINE ABC 
COMMON 
1/BL17/DEL 

2/BL25/E (4,4) , F ( 4,4) ,G < 4 ,4) 

3/BLl/A(4,4) ,BEE(4,4) ,C ( 4 , 4 ) / BL 1 2/TDL I , TDEL 
D2=2. /DEL 
DO 1 1=1,4 
DO 1 J = 1 ,4 
DE I J=D2*fc ( I , J ) 

F I J = F ( I . J) 

BEE( I, J)=-2.*DEIJ+TDEL*G( I, J ) 

C ( I , J) =DEI J-FI J 
1 A ( I , J ) =DEI J+F I J 
RETURN 
END 


SUBROUTINE PANDD ( K, MN ) 

COMMON 

1/IBL4/KMAX,KL 

2/BLl/A(4,4) , BEE (4, 4) ,C(4,4) 

5/BL4/P(4,4,200) ,X(4,200) , ZF1 M ( 4 , 4 , 10 ) , ZF2 M( 4 , 4, 10 ) , ZF3 
4ZF4M (4,4,10) 

1/BL34/DEE(4,4,200) , DST ( 4, 4 , 200 ) 

DIMENSION TM ( 4 , 4 ) , I PIVOT (4) , INDEX (4,2 ) , X2 (4) 
IKL=K+KMAX*(MN-1 ) 

KL I = I KL-1 
DO 1 1=1,4 
DO 1 J=l,4 
SUM=0. 

DO 2 L = 1 ,4 

2 SUM=SUM+C( I ,L)*P(L,J,KLI) 

1 TM(I ,J)=BEE(I,J)-SUM 

CALL MAT INV(TM, 4, X2,0, DETERM, I PIVOT, INDEX, 4, I SCALE) 

DO 5 1=1,4 
DO 5 J=1 ,4 
SUMA=0. 

SUMC=0. 

DO 6 L=1 ,4 
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SUMA=SUMA+TM( I ,L)*A(L, J) 

SUMC=SUMC+TM ( I ,L)*C(L, J) 

P(I,J,IKL)=$UMA 

DEE( I,J,IKL)=TM< I , J ) 

DST(I,J,IKL)=SUMC 

RETURN 

END 


SUBROUTINE H J ( K, MN ) 

COMMON 

0/BL8/R ( 200) ,G AM( 200 ) , OMT ( 200 ) 

2/BL20/DE0MXI 200 ) 

3/BL11/0MXI (200) ,PHEE,T0,T2 
3/BL14/LAM2 * LSD18,LSD1N 

4/BL1 5/NU , Ui ( 10 ) ,V1(10) ,Wl(10),V2( 10) ,U2( 10) ,W2(1G),U3< 
5 W3 ( 10 ) 

6/BL 1T/DEL 

7/BL23/J AY (4,4)»H(4,4) 

8/ IBL2/N ( 10) jMNINIT/ IBL4/KMAX, KL 
EQUIVALENCE ( L2 , LAM2 ) 

REAL L2 , LAM2,LSD1N, LSD18,JAY,NU 

CALL BDB(K,B,0B, D, DD) 

YAH=1 • 

IF(K.EQ.1.0R.K.EQ.KMAX)YAH=2. 

Dl=( l.-NU) 

GA=GAM ( K ) 

OX=OMX T ( K ) 

RA=R ( K) 

EN=N( MN ) 

ENR=EN/RA 

REG=0. 

I F ( YAH.EQ.2. ) REG= 1 • 

OT=OMT ( K ) 

OXT=3.*OMXI (K) -OMT(K) 

0TX=3**0MT(K)-0MXI (K) 

DL=D*L2*D1*ENR 
H ( 1 » 1 ) =B 
H( 1,2)=0. 

H ( 1*3) =0 • 

H ( 1 t A) =0* 

H( 2,1)=0. 

H(2,2)=B*Dl/2.+L2*D*Dl/8.*0TX**2*REG 

H(2,3)=DL/2« *GTX*REG 

H<2,4)=0. 

H(3,l)=0. 

H( 3,2 )=DL*OTX* YAH/4* 

ENR2=ENR**2 

H(3,3)=L2*D*Dl*(YAH*ENR2+( 1 . +NU ) *GA**2 ) 

GA2=GA**2 

H(3,4)=L2 

H(4,l)=0. 

H( 4 , 2 ) =0# 

H<4,3)=-1. 

H(4,4)=0. 

JAY (1,1 )=NU*GA*B 
JAY(1, 2 )=NU*B*ENR 
JAY( 1,3)=B*(GX+NU*0T) 

JAY(1,4)=0. 

JAY (2,1 )=-B*Dl*ENR/2.-DL/8.*OXT*OTX*REG 

J AY ( 2 , 2 ) =-GA*H (2,2) 

JAY(2,3)=-GA*H<2,3) 

J AY( 2 , 4 ) =0* 

JAY (3,1 ) =— L2*D*D1* ( ( 1 • +NU ) *GA2*OX+ENR2/ 4* *OXT*YAH ) 

JAY (3, 2) =-GA*DL/2.*( 2.+0T* ( 1 . +NU ) +0TX/2 .*YAH ) 

J AY ( 3, 3 ) =-L2*D*Dl* ( l.+NU+YAH ) *GA*ENR2 
J AY( 3 , 4 ) =L2*D1 *GA 
JAY (4,1) =0X 
J AY ( 4, 2 ) =0 • 

J AY ( 4, 3 ) =0* 

JAY ( 4,4 ) =0» 

DO 1 1=1,4 
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DO 1 J = l,4 

1 H(I,J)=H(I , J)/2./DEL 
RETURN 
END 


SUBROUTINE TEAETA(K) 

REAL NU,MT 

COMMON/ BL5/TT ( 10) ,MT( 10) ,DT< 10) ,DMT(10) 

1/IBL1/MNMAX 
2/IBL2/N( 10) tMNINIT 

3/IBL7/MNMAX0,MAXD< 10 ) , M AXS ( 10 ) , MAXSY ( 10 ) * I S ( 10 * 10 ) * JS ( 
4*10) *JD< 10*10) T I JS (10 ) 

1/IBL8/LSTEP t ITR 
2 / IBL13/ ITRMAX* LSMAX 
8/BL6/Z( 4,220) , SO E * OSE T ALOAD 
6/BL7/D1 ,S1 

0/BL8/R(200) , GAM (200) ,OMT(200) 

8/BL12/TDLI,TDEL 

9/BL15/NU»Ul( 10 ) ,V1(10) *W1<10) ,V2<10) ,U2<10) ,W2(10),U3< 
0W3 ( 10 ) 

1/BL27/BX3( 10) ,BT3(10) ,BXT3 < 10 ) t BE3 ( 1C ) 

2/BL10/PHIX( 10) , PHIT(IO) ,PHI( 10) 

3/BL28/EXX3( 10 ) * ETT3 (10) * ETX3 ( 10) * EXT3 ( 10) * E X3 ( 1C ) ,ET3( 
3/BL11/0MXK200) ,PHEE,T0,T2 
DIMENSION TX(10) *TTH(10) ,TXT( 10) 

RRA=1 • / R ( K ) 

GA=GAM ( K ) 

OX=OMX I ( K ) 

0 T = 0 MT ( K ) 

CALL BDB(K,BStDB,DS,DD) 

DO 1 M= 1 * MNMAXO 
EN=N( M) 

CALL TLOAD(K) 

TTS=TT ( M ) *ALOAD 

EX=(U3( M)-Ul(M) )*TDLI+OX*W2(M)+OSE*(BX3(M )+BE3(M) ) 

ET = EN *V2(M)*RRA+GA*U2(M)+0T*W2(M)+0SE*(BT3(M)+BE3(M) 
EXT=.5*(TDLI*< V3(M)-V1(M) >- EN *U2 ( M ) *R RA-GA* V2 (M) )+0S 
TX<M) = BS*< EX+NU*ET)-TTS 
TTH(M) = BS 3 M ET + NU* E X ) -TTS 

1 TXT(M)=BS*D1*EXT 
DO 9 M=1 j MNMAX 
SMF=0. 

SMS=0. 

SMV=0. 

SME=0. 

SMN=0* 

SMT=0. 

IF(N(M) •EQ.O) GO TO 20 
MAXL=MAXS ( M ) 

I F ( MAXL • EQ* 0 ) GO TO 2 
DO 3 L=1 *MA XL 
I=IS<L*M) 

J = J S ( L t M ) 

SMF=SMF+TX( I )*PHIX(J)+TX(J)*PHIX( I ) 

SMS=SMS+TTH( I )*PHIT< J)+TTH( J)*PHIT(I ) 

S MV=SMV-PH I T ( I ) *TXT ( J ) -PH IT ( J ) *TXT ( I ) 

SME=SME+PHI X ( I )*TXT(J )+PHIX< J)*TXT( I ) 

S MN=SMN + TX ( I )*PHI(J)+TX(J)*PHI ( I ) 

3 SMT=SMT+TTH(I)*PHI ( J ) +TTH ( J ) * PH I ( I ) 

2 MAXL=MAXD( M ) 

IF(MAXL.EQ.O) GO TO 4 
00 5 L= 1 * MAXL 
I-ID(LfMi 

J=JD( L,M) 

$MF=SMF+TX ( I )*PHIX( J)+TX< J)*PHIX< I ) 

SMS=SMS— TTH{ I ) *PH IT (J ) + TTH(J) *PH IT ( I ) 

SMV=SMV+PHIT (I )*TXT(J)+PHIT (J )*TXT ( I ) 

S ME=$ ME-PH I X ( I )*TXT(J)+PHIX(J)*TXT( I ) 

S MN=SMM-TX ( I )#PHI(J)+TX(J)*PHI ( I ) 

5 SMT=SMT-TTH( I )*PHI ( J ) +TTH ( J ) * PH I (I) 

4 IF< MAXSY(M) .EQ.O ) GO TO 10 
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I = 1 JS C M ) 

SMF=SMF+TX< I)*PHIX(I1 
SMS=SMS+TTH( I ) *PHIT ( I ) 
SMV=SMV— PHIT ( I ) *TXT ( I > 
SME=SME+PHIX( I )*TXT(I ) 
SMN=SMN+TX ( I )*PHI ( I ) 
SMT=SMT+TTH(I)*PHI (I) 

GO TO 10 

20 DO 21 L=1 , MNMAXO 

„ SMF=SMF+TX(L)*PHIX<L) 

21 SMV=SMV+PHIT(L )*TXT(L) 
IFIM.GT. MNMAXO) GO TO 10 

„ SMF=SMF+TX( M)*PHIX( M) 

10 EXX3(M)=SMF*»5 
ETT3(M)=SMS*.5 
ETX3(M)=SMV*.5 
EXT3< M)=SME*.5 
EX3 ( M) = SMN* • 5 
ET3(M)=SMT*.5 
9 CONTINUE 

DO 30 M=l, MNMAXO 
U1(M)=U2(M» 

VI < M) =V2 ( M ) 

W1(M)=W2(M) 

U2(M)=U3(M) 

V2(M)=V3(M) 

30 W2(M)=W3(M) 

RETURN 

END 


SUBROUTINE FORCE ( K ) 

REAL NU,MT,LAM2, MASS, MAS 
COMMON 

3/IBL4/KMAX,KL 

i^it!^ L ?;^i};E^? L x i2/K " #xiiK " 4x2 ’ NcoNv 

£7lkSf? t ? , ?A? 0 ° , ’ XU ' 200,,ZF1M{4 » 4 * iol * 2F2M(4 ' 4 '10>’ ZF 3 

5 ZF4M( 4, 4, 10 ) 

8/BL6/Z 4,220) , SOE ,OSE , ALOAD 
7/BL7/D1 , SI 

0/BL8/R ( 200 ) , GAM(200),0MT(200) 

9/BL9/FF S ( 4 * 10 ), ELIS (4) ,GEES(4, 10) 

3/BL11/QMXK200) ,PHEE,T0,T2 
1/BL12/TDLI ,TDEL 
2/BL14/LAM2,LSD18,LSDlN 

?^IV}^ NU,ui,10) * V1( 10 > * Wit 10) ,V2( 10) ,U2( 10) ,W2(10),U3( 
4W3( 10) 

5/BL17/DEL 

BL24h/ D L (4,4,10),DG(4,4,10),DF(4,4,10) 

COMMON 

1/BL27/BX3( 10) , BT3( 10) ,BXT3( 10) ,BE3( 10) 

$;iv^fs)U5?iifj{!a?iiS{;m?jjf xiao, ' ETi,io '' Exx2 ' 
S:!? I SJ WWl'l* M § I " 4X1 10: ’■ • PTl 11 101 1 

l/BL102/oIho°D/iLloj;HASsf)o8| ,20 '‘ i ’ 220 '’ Z2 ' 4 ’ 220l,Z3l ‘' 

DIMENSION GEEt 4 ) 

FDIFFl A,B,C)=(-1.5*A+2.*B-.5*C) /DEL 
RS=R ( K ) 

RR=1./RS 
GA=GAM( K ) 

OX=OMXI ( K ) 

OT=OMT ( K ) 

DL2=D1*L AM2 

CALL BDB(K,BS,DBS,D,OD) 
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CALL PLOAD(K) 

CALL TLOAD(K) 

M AS=MAS S ( K ) 

DO 4 M=1,MNMAX 
I Z=K+1 + ( M— 1 )*KMAX2 
I K=K+ ( M—l ) ^KMAX 
I Kl= I K- 1 
EN=N( M ) 

ENR=EN*RR 
EMT = MT ( M ) 

GEE ( 1 ) = (-PX( M)+DT(M)-DL2*GA*OX*EMT)*TDEL*ALOAD 

1 + MA$*<-5.*ZO<l,IZ )+4.*Z2(l, IZ)-Z3<1,IZ) )*TDEL/DELS 

GEE( 2 ) = <-PT< M)-ENR*TT(M)-DL2*ENR*OT*EMT)*TDEL*ALOAD 

1 + MA$*(-5»*ZO<2, IZ ) +4* *Z2 ( 2 T IZ)-Z3(?*IZ) ) *TDEL /DELS 

GEE ( 3 ) = ( — PR < M)-< OX+OT)*TT< M)-DL2* ( GA*D MT < M ) - ( OX*OT- 

1*MT(M) ) )*TDEL*ALOAD 

1 + MAS*(-5.*ZO(3, IZ ) +4. *Z2 ( 3 * I Z )-Z3 ( 3 * IZ) )*T DEL/DELS 

l GEE ( 4 ) = MT(M)*TDEL*ALOAD 
I F ( ITRMAX.EQ.l ) GO TO 50 
IF(K.GT.l) GO TO 6 
BX2T=BX1 ( M ) 

BT2T=BT1 ( M ) 

BXT2T=B XT1 ( M) 

BE2T=BE1 ( M) 

EXX2T=EXX1 ( M ) 

ETT2T=ETT1 ( M ) 

ETX2T=ETX1(M) 

EXT2 T=E XT1 ( M ) 

EX2T=EX 1(M) 

ET2T=ET1 { M ) 

DBX= FDIFF(BX2T,BX2(M) ,BX3(M) J 
DBT = FDIFF(BT2T,BT2(M) ,BT3(M) ) 

DBXT= FDIFF(BXT2T,BXT2(M) ,BXT3(M) ) 

DBE= FDIFF(BE2T,BE2(M) ,BE3(M) ) 

DET = FDIFF(ET2T,ET2(M) , ET3 ( M ) ) 

DEX= FDIFF(EX2T,EX2(M) ,EX3(M) ) 

DEXX= FDIFF(EXX2T,EXX2(M) ,EXX3(M) ) 

DETX= FDIFF(ETX2T, ETX2(M) f ETX3(M) ) 

GO TO 7 

6 IF(K.LT.KMAX) GO TO 8 
BX2T=BX3(M) 

BT2T=BT3 < M) 

BXT2T=B XT3 < M ) 

BE2T=BE3 ( M ) 

E XX2T=E XX3 ( M ) 

ETT2T=ETT3(M) 

ETX2T=ETX3(M) 

E XT2T=E XT3 ( M ) 

EX2T=EX3 ( M ) 

ET2T=ET3 ( M) 

DBX=-FDIFF(BX2T,BX2(M) ,BX1(M) ) 

DBT=-FDIFF(BT2T » BT2 ( M ) * BT1 ( M ) ) 

DBXT=-FDIFF (BXT2T, BXT2 (M) * BXT1 (M) ) 

DBE=-FD I FF ( BE2T t BE 2 ( M ) , BE1 ( M ) ) 

DET=-FDIFF( ET2T,ET2(M) t ET1(M) ) 

DEX=-FDIFF(EX2T t EX2(M) ,EX1(M) ) 
DEXX=-FDlFF(EXX2T t EXX2(M) , EXX1(M) ) 
DETX=-FDIFF(ETX2T,ETX2(M) , ETX1(M) ) 

GO TO 7 

8 DBX=TDL I# ( BX3< M)-BXl(M) ) 

DBT=TDL I*< BT3( M) -BT1 ( M ) ) 

DBE=TDL I*(BE3(M)-BE1(M) ) 

DBXT=TDLI*( BXT3 (M)-BXTl (M) ) 

DET=TDLI*(ET3( M)-ETl(M) ) 

DEX=TDLI*<EX3( M)-EXl(M) ) 

DEXX=TDLI*( EXX3 (M)-EXXl (M) ) 

DETX=TDLI*< ETX3(M)-ETX1(M) ) 

BX2T=BX2(M) 

BT2T=BT2 ( M) 

BXT2T=B XT2 ( M ) 

BE2T=BE2 { M ) 

EXX2T=EXX2(M) 
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§1 J?I = ETT2 ( M ) 

EXtIt^E XT? (M l 

TjGEEaflc— •’ 

2 

i GEE( 

»Qpf® 

SUMxio.’ M, ” GEe( 11 

DO 21 j-i 4 

Sllfrf * elis,j, ' os "’-— 

SUMxio. ,4 

DO 12 J=1 a 

^ l||llP^^ e( r,J ' IK, * GEE(J, ' DST ' I * J ’ IK, * X( 

arar iNE ^date 

i^8Ll/MNMAX 

israH?i«p»« »» ”r “ 

BMRX2(M) X 
glr j ^> 7 BT 2 ( M) 

ffifflifliwr 

!||<«!: § r?!S! 
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IrrifM 

ETX1 8 iffjf JJJ 
fXTl(M =ExrIrM 

flfl ( M ) =IttI ( M ) 

||2 T fAri|| 3 x ^r» 

RE?A^ =ET3 <«» 
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SUBROUTINE MAT I NV ( A , N , B , M * DETERM, I P I VOT , I NDE X, NMAX, I SC 
C MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR 

C 
C 

DIMENSION I PI VOT IN) , A ( NMAX , N ) , B < NMAX , M ) , I NDE X ( NMA X , 2 ) 
EQUIVALENCE (IROW,JROW), ( ICOLUM t JCOLUM ) , (AMAX, T, SW 

C INITIALIZATION 

C 

5 I SCALE=0 

6 R 1= 10. 0**1 8 

7 R2=1.0/R1 
10 DETERM=1.0 
15 DO 20 J = 1 , N 
20 I PI VOT ( J )=0 
30 DO 550 I =1 , N 

C 

C SEARCH FOR PIVOT ELEMENT 

C 

40 AMAX=0.0 
45 DO 105 J = 1 , N 

50 IF (IPI VOT ( J ) - 1 ) 60, 105, 60 
60 DO 100 K=1 , N 

70 IF ( IPI V0T(K)-1 ) 80, 100, 740 

80 IF (ABS(AMAX)-ABS(A( J,K) J >85,100,100 

85 IROW=J 

90 I COLUM=K 

95 AMAX=A ( J , K ) 

100 CONTINUE 
105 CONTINUE 

110 I PIVOT ( ICOLUM) = IPI VOT ( ICOLUM )+l 
C 

C INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

C 

130 IF ( IROW-ICOLUM) 140, 260, 140 
140 DETERM=-DETERM 
150 DO 200 L=ltN 
160 SWAP=A ( I ROW , L ) 

170 A( IROW,L)=A( ICOLUM, U 
200 A ( ICOLUM, L)=SWAP 
205 I F ( M ) 260, 260, 210 
210 DO 250 L=1 , M 
220 SWAP=B ( I ROW , L ) 

230 B ( I ROW ,L)=B( ICOLUM, L) 

2 50 B( ICOLUM, L) = SWAP 
260 I NDEX { 1 , 1 )= I ROW 
270 I NDEX ( 1,2) = ICOLUM 
310 PIVOT=A( ICOLUM, ICOLUM) 

C 

C SCALE THE DETERMINANT 

C 

1000 P I VOT 1 = P I VOT 

1005 I F ( ABS ( DETERM ) -R 1 >1030,1010,1010 
1010 DETERM=DETERM/R1 
ISCALE=ISCALE+1 

IF( ABS ( DETERM )-Rl) 1060, 1020,1020 
1020 DETER M= DETERM/ R 1 
ISCALE=ISCALE+1 
GO TO 1060 

1030 IF( ABS( DETERM) -R2) 1040, 1040,1060 
1040 DETERM=DETERM*R1 
ISCALE=ISCALE-1 

IF( ABS( OETERM)-R 2) 1050, 1050, 1060 
1050 DETERM=DETERM*R1 
ISCALE= ISCALE-1 

10 60 I F ( AB S ( PI VOT I ) —R1 ) 1 0 90 , 1070 , 1070 
1070 P I VOT I =P I VOTI / R 1 
ISCALE=ISCALE+1 
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non non non 


IF( ABS< PIV0TIJ-R1 >3 20,1080, 1080 
1080 P IVOT I=PIV0TI/R1 
I SCALE= I SCALE + 1 
GO TO 320 

1090 I F ( ABS( PIVOT I )-R2) 2000, 2000, 3 20 
2000 PIV0TI=PIV0TI*R1 
I SC ALE= I SCALE- 1 

IF( ABS( PIVOT I )-R2 >2010,2010, 320 
2010 P I VOT I=P I VOTI *R1 
ISCALE=ISCALE-1 
320 DETERM=DETERM*PIVOTI 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 A ( ICOLUM , I COLUM ) =1 • 0 
340 DO 350 L = 1,N 

350 A ( ICOLUM, L)=A( ICOLUM, L ) /PI VOT 
355 I F ( M ) 380, 380, 360 
360 DO 370 L=1 , M 

370 B( ICOLUM, L)=B( ICOLUM, L ) /PIVOT 

REDUCE NON-PIVOT ROWS 

380 DO 550 Ll=l , N 

390 I F ( LI— I COLUM ) 400, 550, 400 

400 T=A( LI, ICOLUM) 

420 A(L1, ICOLUM )=0*0 
430 DO 450 L=1,N 

450 A(L1,L)=A(L1,L)-A( IC0LUM,L)*T 
455 I F ( M ) 550, 550, 460 
460 DO 500 L=1 , M 

500 B(L1,L)=B(L1,L)-B( ICOLUM, L)*T 
550 CONTINUE 

INTERCHANGE COLUMNS 

600 DO 710 1=1, N 
610 L=N+ 1— I 

620 IF ( INDEX(L,1)-INDEX(L,2> ) 630, 710, 630 
630 JROW= INDEX(L,1) 

640 JCOLUM=INDEX(L ,2 ) 

650 DO 705 K=1 , N 
660 SWAP = A( K, JROW ) 

670 A(K,JROW)=A(K, JCOLUM) 

700 A(K,JCOLUM)=SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 


SUBROUTINE PMATRX 

REAL JAY 

COMMON 

1 / I BL1/MNMAX 

2/IBL2/N(10) , MN I NIT 

3/IBL3/M0,Ml,M2,M3 

4/IBL4/KMAX,KL 

5/IBL5/IBCINL,IBCFNL 

6/BLl/A( 4,4),BEE(4»4),C(4,4) 

5/BL4/P(4,4,200) , X( 4,200 ) , ZF1M( 4, 4, 10 ) , ZF2 M( 4, 4, 10), ZF3 
8 ZF4M ( 4,4 , 10 ) 

9/BL13/0MEGl(4,4) ,CAPL1(4,4) , OMEGL ( 4 , 4 ) , C A PLL ( 4, 4 ) ,UNIT 
A/BL23/JAY(4,4) ,H(4,4) 

B/BL24/DL( 4,4,10) , DG ( 4 , 4 , 10 ) , DF ( 4 , 4 , 10 ) 

C/BL25/E(4,4) ,F(4,4) ,G(4,4> 

DIMENSION PAT A ( 4 , 4 ) ,PBTA(4,4) , POT A (4, 4) , PJTA<4,4) , DLL ( 
14) ,DGG< 4,4) ,ZF 1(4,4 ),ZF2( 4,4) ,ZFPO( 4,4) ,ZFP1(4,4) ,ZFP2 
2T(4), INDEX (4,2 ) , CLO ( 4 , 4 ) , CL 1 ( 4, 4) ,CL2(4,4) ,G1 (4) 
EQUIVALENCE (CLOU ) ,ZF1M(1) ) , (CL1(1 ) , ZF2M ( 1 ) ) »(CL2( 1), 
KZFPO(l), PAT A ( 1 ) ), (ZFP1(1),PBTA( 1) ) ,(ZFP2(1) ,P0TA(1) ) 
2, (ZF1(1) , DLL (1) ), ( ZF2 ( 1 ) , PTR ( 1 ) ) 
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IF< IBCINL.LT.O) GO TO 10 
DO 1 MN=MNINIT , MNMAX 
CALL HJ(1,MN) 

CALL EFG(1,MN) 

CALL ABC 

CALL MATINV < C,4,Gl,0t DETERM f I PIVOT, INDEX ,4, ISCALE ) 

DO 3 1 = 1,4 
DO 3 J= 1 1 4 
SUMA=0* 

SUMB=0# 

SUMO=0. 

SUMJ=0. 

DO 4 L= 1 , 4 

SUMJ=SUMJ+0MEG1 <ItL)*JAY(L, J) 

SUMA=SUMA+C (I , L ) *A ( L, J ) 

$UMB=SUMB+C (I.L)*BEE(L,J) 

4 SUM0=SUM0+0MEG1(IjL)*H(LtJ) 

PATA(I , J)=SUMA+UNIT<I , J) 

PBTA(I, J)=SUMB 

POT A ( I , J )=SUMO 
3 P JTA ( I , J ) = SUMJ 
DO 5 1=1,4 
DO 5 J=1 ,4 
SUM0B=0. 

SUMOA=0. 

SUMOC=C. 

DO 6 L= 1 1 4 . 

SUMOB=SUMOB+POTA(I , L)*PBTA(L, J) 

SUMOA=SUMOA+POTA( I , L ) *PATA ( L , J ) 

6 SUMOC=SUMOC+POTA(I, L)*C<L,J) 

DLL ( I, J)=SUMOB + PJTA( It JJ+CAPLK I, J) 

PTR ( I j J ) =SUMOA 

5 DGG( I, J)=SUMOC 

CALL MATINV(DLL,4,PTR,4,DETERM, I P I VOT 1 1 NDE X , 4 1 1 SC AL E ) 
DO 1 1 = 1 1 4 
DO 1 J = 1 ,4 
SUMD=0* 

SUME=0« 

DO 7 L = 1 1 4 

SUMD=SUMD+DLL( I , L ) *DGG < L , J > 

7 SUME=$UME-DLL( I , L ) *0MEG1 < L , J ) 

DL(I,J,MN)=DLL<I,J) 

DG(I, J»MN)=SUMD 
DF< I, J,MN)=$UME 
I J=1+KMAX*( MN-1 ) 

1 P ( I ,J,I J)=PTR( I , J) 

GO TO 20 

10 DO 11 MN=MNINIT, MNMAX 
NN=N( MN ) 

IJ=1+KMAX*<MN-1 ) 

DO 14 1=1,4 
X(ItIJ)=0, 

DO 14 J=l,4 
14 PUtJ,IJ>=0. 

I F ( NN «GT • 3 ) GO TO 11 
IF(NN.GT.2> GO TO 90 
IF(NN.GT.l) GO TO 12 
IF(NN.GT.O) GO TO 13 
P ( 3 1 3 , I J )=-l • 

P(4,4,I J)=-l. 

MO=MN 
GO TO 11 

12 P(4,4tl J)=-l. 

M2=MN 

GO TO 11 
90 M3=MN 
GO TO 11 

13 PU,1,IJ)=-1. 

P ( 2, 1 , IJ ) = 1 • 

M1 = MN 

11 CONTINUE 
20 KLAST=KMAX 
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23 


42 

41 


44 

43 


46 

45 

40 

30 


35 


IF( IBCFNL.LT. 0) KL AST=KL 
DO 23 K=2, KLAST 
DO 23 MN=MNINI T , MNMAX 
CALL EFG(K,MN) 

CALL ABC 

CALL PANDD(K,MN) 

IF( IBCFNL.LT. 0) GO TO 30 
DO 40 MN=MNINIT, MNMAX 
IKL=MN*KMAX-1 
JKL=KMAX*MN 
CALL HJ { KMAX, MN ) 

DO 41 1=1,4 
DO 41 J = l,4 
SUM0=0. 

$UMP=0. 

SUM J=0 • 

DO 42 L = 1 , 4 

SUMO= SUMO+OMEG L ( I , L ) *H ( L , J ) 

SUMP=SUMP+P< I.L, IKL)*P(L, J, JKL) 
SUMJ=SUMJ+OMEGL (I , L)*JAY(L, J ) 

PATAU , J) = SUMO 

PBTAI I, J )=UNIT ( I t J ) -SUMP 

PJTA(I t J)=SUMJ+CAPLL(I , J) 

DO 43 1=1 t 4 
DO 43 J = 1 1 4 
$UM0P=0. 

SUM J P=0 • 

SUMOM=0 • 

DO 44 L = 1 1 4 

SUMOP=SUMOP+PATA<I ,L) *PBTA(L, J) 
SUMJP=SUMJP+PJTA(I ,L) *P<L,J, JKL) 
SUMOM=$UMOM-PATA(I , L)*P(L f J,IKL> 

ZF1 ( I t J ) =$UMOP-SUM J P 
ZF2( I , J)=SUMOM-PJTA<I , J) 

CALL MATINV( ZF 1 , 4 t ZF2 ,4,DETERM, I PI V0 T t I NDEX,4 
DO 45 I = 1 1 4 
DO 45 J = 1 1 4 
S ZF3=0 • 

S ZF4=0 • 

DO 46 L = 1 1 4 

SZF3 = SZF3 + ZF1( I , L )*PATA(L, J ) 

SZF4=SZF4-ZFl ( I t L ) * OME GL ( L , J ) 

ZF3M ( I , J ,MN)=SZF3 
ZF4M ( I , J,MN)=SZF4 
ZF1 M ( It JtMN)=ZFl (I T J> 

ZF2M ( I, J,MN)=ZF2(I , J) 

CONTINUE 

RETURN 

DO 31 MN=MNINIT, MNMAX 
I KL= MN*KMAX-1 
NN=N( MN ) 

IF ( NN .GT.3) GO TO 31 
I F ( NN.GT .2 ) GO TO 300 
I F ( NN .GT.l) GO TO 33 
I F ( NN .GT.O) GO TO 34 

MO = MN 

DO 35 J = 1 1 4 
DO 35 1=1,4 
CLO ( I, J)=0. 

ZFPO( I , J)=0. 

ZFP0(1,1)=1. 

ZFP0(2,2)=1. 

ZFPO( 3,1)=P(3,1,IKL) 

ZFP0(3,2)=P(3,2, I KL ) 

ZFP0(3,3)=P(3t3, IKL )+l • 

ZFPO( 3,4)=P( 3,4, IKL) 

ZFP0(4, 1 )=P(4, 1, IKL) 

ZFP0(4,2)=P(4,2 , IKL) 

ZFP0(4,3)=P(4,3, IKL) 

ZFPO (4,4)=P(4,4, IKL )+l. 

CL0( 3 , 3) =1. 

CLO ( 4 , 4 ) = 1 • 
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300 

34 


60 


33 


70 


31 


CALL MATINV<ZFP0,4,CL0, 4, DETERM, I PIVOT, INDEX, 4, 1 SCALE ) 
GO TO 31 


M3=MN 
GO TO 
M1=MN 
DO 60 
DO 60 


31 


J = 1 , 4 
__ 1=1,4 
CL1 ( I , J ) =0* 

ZFP1( I, J )=0. 

ZFP1(1,1)=P(1,1,IKL)+1. 

ZFP1(1,2)=P(1,2, IKL) 

ZFPl(l,3) = P(l,3,T KL ) 

ZFP1(1,4)=P(1,4,IKL) 

ZFPK2, 1) = 1. 

ZF PI ( 2 , 2 )=-l • 

ZFP1 ( 3 , 3 ) -1 • 

Z Li(i 4 it-i 1 * 

CAL^AATINVJZFPl^fCLl ,4 , DETERM, I PI VOT, I NDEX , 4, 1 SCAL E ) 

GO TO 31 

M2=MN 

DO 70 J = 1 , 4 
DO 70 1=1,4 
CL2 ( I • J ) =0. 

ZFP2 ( I , J >=0. 

ZFP2 ( 1 , 1 ) = 1 • 

ZFP2( 2 , 2 )=1. 

ZFP2( 3 , 3 )= 1. 

ZFP2(4,1)=P(4,1,IKL) 

ZFP2(4,2)=P(4,2,IKL) 

ZFP2(4,3)=P(4,3, IKL) 

ZFP2(4,4)=P(4,4, IKL1+1. 

0 L 2 ( j — 

CALL MAT I NV (ZFP2*4,CL2 , 4, DETERM , I PI VO T, INDEX, 4, 1 SCALE 1 

CONTINUE 

RETURN 

END 


SUBROUTINE POLE(K) 

REAL NU , MT , MX, MTH, MXT , MTS , KX , KT , KXT , LAM , L AM2 , M AS S 
COMMON /IBL2/NU0) ,MNINIT 
2/IBL3/M0,Ml,M2,M3 
1 / 1 BL4/KMAX, KL 
2/IBL5/IBCINL,IBCFNL 

3/IBL7/MNMAX0,MAXD( 10) ,MAXS( 10) , MAXSY ( 10 ) , IS ( 10 ,1C » , JS( 
4,10) , JDC 10, 10) , I JS( 10 ) 

5/IBL8/LSTEP,ITR 
6/IBL10/I FREQ, NTH MAX 
7/IBL12/KMAXl,KMAX2,NC0NV 
8/BL6/Z<4,220) , SOE, OSE , ALOAD 
9/BL7/D1 , SI 

0/BL8/R(200) ,GAM(200),QMT(20C) 

3/BL11/0MXK200) ,PHEE,T0,T2 
4/BL20/DEOMX(200) 

1/BL10/PHIX( 10) , PHI T ( 10 ) , PHI ( 10 ) 

3/BL12/TDLI »TDEL 
4/BL14/LAM2,LSD18,LSDlN 

5/BL15/NU ,U1 (10) , VI ( 10) ,W1( 10) ,V2( 10) ,U2( 10) ,W2( 1C),U3( 
6W3( 10) 

8/BL27/BX3(10)»BT3( 10 ) , B XT3( 10 ) , BE 3( 10 ) 

COMMON/ BL32/TKN, EL AST, CHAR, SI GO 
1/BL31/DELSQ,EXT1(10) 

2/BL17/DEL 
3/BL19/TH ( 6 ) 

COMMON/ BL5/TT( 10 ) , MT( 10 ) , DT( 10 ) , DMT ( 10) 

COMMON/ BL4/ P( 4,4,200) , X ( 4 . 200 ) , ZF 1M< 4,4, in ) , ZF2M( 4, 4, 1 
1»ZF3M(4,4,1C) ,ZF4M(4,4,10) 

COMMON /BL100/S0RD,TFE0/BL101/Z0( 4,220) ,Z2(4,220) ,Z3(4 
1/BL102/DEL0AD/BL103 /MASS (200) 

2/BL110/TX( 10) ,TTH( 10) ,TXT( 10) ,MX( 10 ) , MTH ( 10) ,MXT( 10 ) ,Q 
3/BL111/ABZ, ABZO, ABZN, ABZ3,DD2 
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CALL BDB(K,BS,DB,DS,DD) 

IF(K.EQ.KMAX) GO TO 301 
DO 202 MN=1,MNMAX0 
U1(MN)=U2(MN) 

VI ( MN) =V2 ( MN ) 

W1(MN)=W2(MN) 

1 3=3+( MM— 1 ) *KMAX2 
12=13-1 

U2(MN)=Z(1, 13) 

V2 ( MN )=Z(2»I3) 

W2(MN)=Z(3, 13) 

PHI X ( MN ) =0. 

PHI T(MN)=0. 

PHI<MN)=0. 

MX(MN)=Z(4,I2)*ABZ3 
MTH(MN) =0. 

MXT ( MN ) =0 • 

QS(MN)=0. 

TX(MN)=0. 

TTH( MN ) =0. 

202 TXT ( MN ) =0« 

IF(Ml.EQ.O) GO TO 203 
CALL INLPOL 

PHIX(M1 ) =PHEE* ABZO 
PHIT (Ml )=-PHEE*ABZO 
1 1=3 + ( Ml— 1 ) =*KM AX2 
IP1=II+1 
I Ml = 1 1- 1 
GA2=GAM( 2) 

CALL BDB(^2jBS,DB,DS,DD) 

PHIXX=-Z(3, IP1 >*TDLI+OMXI (2)*Z< 1, 1 1 ) 
PHITT=Z(3,II)/R(2)+0MT(2)*Z(2,I!) ,, IT , ,,, 

PHII=((Z(2,IP1)-Z(2,IM1) ) *TDL I+GA2*Z ( 2, II)+Z(1,II)/R<2 
QS(M1)=SIG0*TKN*LAM2*( ( 2 .-NU ) *Z ( 4 ,1 1 )-DS*DD2* ( PHI TT/R ( 

1 *GA2 )+Dl*MT< Ml ) * ALOAD+DS*D 1* ( -PH I XX/R ( 2 ) -GA2*PHI TT+ 

2 (2))*PHII + (Z( 3, I PI )*TDLI-GA2*Z< 3, II ) )/R<2)+GA2*< 

3 < 2) )*Z<2, I I ) +OMT (2)*(Z(2,IP1)-Z(2»IM1) )*TOLI )*.5) 
IF(MO.EQ.O) GO TO 204 

TX(MO) =TC* AB Z 
TTH(MO) =T0#AB Z 
MTH(MO) =MX ( MO ) 

204 IF(M2.EQ.O) GO TO 205 

TX ( M2 ) = T 2*AB Z 
TTH( M2 ) =-T2*ABZ 
TXT(M2)=-T2*ABZ 
MTH( M2 ) =-MX ( M2 ) 

MXT ( M2 ) =MTH(M2 ) 

GO TO 205 

203 IF(MO.EQ.O) GO TO 206 
I3=3+( MO-1 ) *KMAX2 
14=13+1 

TX^MoI=BS*Sl* ( <2.*Z<1,I3)-.5*Z<1, 14) ) /DEL+OMX I ( 1 ) *Z < 3 , 
1 -TT ( MO) *ABZ*ALOAD 

TTH( MO ) =TX ( MC' ) 

MTH( MO) =MX( MO) 

206 I F ( M2 • EQ.O ) GO TO 205 
I3=3+(M2-1)*KMAX2 

14=13+1 

TX(M2)=BS*D1*(2.*Z( 1, 13 )-.5*Z( 1,14) )/DEL 
TX( M2 ) = TX(M2).*ABZ 
TTH( M2 ) =— T X ( M 2 ) 

TXT ( M2 ) =— TX ( M2 ) 

MTH ( M2 ) =— MX ( M2 ) 

MXT ( M2 ) =-MX( M2 ) 

205 RETURN 
301 CONTINUE 

DO 302 MN=1 , MNMAXO 
U1 ( MN ) =U2 ( MN ) 

VI ( MN ) =V2 ( MN J 
W1(MN)=W2( MN ) 
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PH IX ( MN ) =0* 

PHIT ( MN ) =0# 

PHI(MN)=0. 

IK=KMAX1+(MN-I )*KMAX2 
MX(MN)=Z(4,IK)*ABZ3 
MTH ( MN ) =0* 

MXT < MN ) =0« 

QS ( MN ) =0« 

TX ( MN ) =0 • 

TTH ( MN ) =0. 

302 TXT(MN)=0. 

IF(Ml.EQ.O) GO TO 303 
CALL FNLPOL 

PHIX1M1 )=PHEE*ABZO 

PH I T ( Ml )= PHEE*ABZO 

I I=KMAX+(M1-1)*KMAX2 

IP1=II+1 

IM1=II-1 

GAK^GAM ( KL ) 

CALL BDB(KL,BS,DB f DS,DD) 

CALL TLOAD(KL) 

PHIXX = Z(3,IM1)*TQLI+0MXI (KL)*Z(1» I I ) 
PHITT=Z(3tII)/R(KL)+0MT(KL)*Z(2,I I ) 

PHII = ( < ZC2.IP1 )-Z(2 t IMl n*TDLI+GAK*Z(2T I I ) +Z ( 1 » 1 1 ) /R( K 
QS(M1)=-SIG0*TKN*LAM2*( < 2.-NU) *Z (4, 1 1 ) -DS*DD2* ( PH ITT /R 
1*GAK ) +D1 *MT ( Ml )*ALOAD-DS*Dl s M-PHIXX/R(KL)-GAK*PHITT+(0 

2 (KL) )*PHII+(-Z<3t IM1)*TDLI-GAK*Z(3,II ) )/R(KL)+GA 

3 -OMT(KL) )*Z(2 f II )+OMT(KL)*< Z(2 , I PI )-Z ( 2 , IM1 ) ) 

4 /DEL 
IF(MG.EQ.O) GO TO 304 

TX(MO) =TO*ABZ 
TTH(MO) =TO*ABZ 
MTH ( MO) =M X ( MO ) 

304 IF(M2.EQ.O) GO TO 305 

TX ( M2 ) = T2*A BZ 
TTH( M2 ) =-T2*ABZ 
TXT ( M2 ) = T2*ABZ 
MTH( M2 ) =— MX ( M2 ) 

MXT ( M2 ) =MX ( M2 ) 

GO TO 305 

303 IF(MO.EQ.O) GO TO 306 
I KM=KMAX+ ( MO-1 ) *KMAX2 
I Ml = I KM-1 

CALL TLOAD(KMAX) 

TX(M0)=BS*S1*( <-2.*Z( 1, IKM)+.5*Z( It IM1) )/DEL+QMXI ( KMAX 
1 )*ABZ-TT(MO)*ABZ*ALOAO 

TTH(MO) =TX ( MO ) 

MTH(MO) =MX ( MO ) 

3C6 I F ( M2 • EQ# 0 ) GO TO 305 
IKM=KMAX+(M2-1)*KMAX2 
I Ml= I KM— 1 

TX(M2) = BS*Dl : M-2#*Z(ltIKM)+* 5*Z ( 1 , 1 Ml ) ) /DEL 
TX ( M2 )=TX(M2)*ABZ 
TTH( M2 ) =— T X ( M2 ) 

T XT ( M2 ) = TX ( M2 ) 

MTH( M2 ) =— MX ( M2 ) 

MXT ( M2 ) = MX ( M2 ) 

305 RETURN 
END 


SUBROUTINE PHIBET(K) 

REAL NU 
COMMON 

1/IBL1/MNMAX 
2/IBL2/N( 10) tMNINIT 
3/IBL4/KMAX,KL 

3/IBL7/MNMAX0,MAXD< 1C) ,MAXS( 10) *MAXSY< 10 ) , I S ( 1 0 i 1C ) , JS ( 
4 1 10 ) t JD ( 10 » 10 ) ,IJS(10) 

6/IBL12/KMAXl ? KMAX2 tNCONV 
2/IBL13/ ITRMAX, LSMAX 
8 /BL6/Z ( 4 »220 ) , $0E , OSE f ALOAD 
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0/BL8/R ( 200 ) iGAM( 200 ) , OMT ( 200 ) 

9/BL10/PHIX( 1C ) ,PHIT(10) ,PHI( 10 ) 

3/BL11/0MXI (200) ,PHEE,T0,T2 
1/BL12/TDLI ,TDEL 

2/BL15/NU,Ul< 10 ) t Vl( 10) ,W1( 10) t V2( 10) ,U2( 10) ,W2(10) ,U3( 
3 W3 ( 10 ) 

4/ BL27 / BX3 ( 10),BT3( 10) ,BXT3( 10) ,BE3( 10) 

OX=OMX T ( K ) 

OT -OMT ( K ) 

RRA=1./R(K) 

GA=GAM( K ) 

KP2=K+2 

DO 1 M= 1 » MNM AXO 
EN=N( M ) 

IK=KP2+(M-1) *K MAX2 
U3(M)=Z< It IK) 

V3 ( M) = Z ( 2 » I K ) 

W3 ( M ) =Z ( 3 IK) 

PHIX(M)=-fDLI*(W3(M)-Wl(M) )+0X*U2(M) 

PH I T ( M ) = EN*W2 ( M ) *RR A+V 2 ( M ) *0T 

1 PHI (M)=(TDLI*(V3( M)-Vl(M) ) +G A*V2 ( M ) +EN*U2 ( M ) * RRA ) *. 5 
IF( ITRMAX.EQ.l ) RETURN 

DO 9 M=1 tMNMAX 
SM0=0 • 

SMT=0« 

SMR=0. 

SMF=0* 

I F ( N< M ) .EQ.O) GO TO 20 
MAXL=MAXS ( M ) 

IF(MAXL.EQ.O) GO TO 2 
DO 3 L= 1 » MA XL 
I=IS(L t M) 

J =JS ( L * M ) 

SMO=SMO+PHI X ( I )*PHIX(J) 

SMT=SMT-PHIT( I )*PHIT( J) 

S MR=S MR+ PH I X ( I)*PHIT( J)+PHIX( J )*PHIT< I ) 

3 SMF=SMF-PHI ( I ) *PHI (J) 

2 MAXL=MA XD( M ) 

I F ( MAX L • EQ • 0 ) GO TO 4 
DO 5 L=1,MAXL 
I =T D ( L * M ) 

J = JD { L ? M ) 

S MO = SMO+ PH I X ( I ) * PH I X < J ) 

SMT=SMT +PH I T ( I )*PHIT( J ) 

S MR = S MR- PH I X ( I )*PHIT( J)+PHIX( J)*PHIT( I) 

5 SMF= SM F+PH I ( I ) * PHI ( J ) 

4 I F( MAXSY(M) .EQ.O ) GO TO 10 
I = I J S < M ) 

SMO=SMO+PHI X ( I )**2/2. 

SMT=SMT-PH IT ( I )**2/2. 

SMR= ( S MR+ PH I X ( I ) *PH IT ( I ) ) 

S MF=SMF- PH I ( I ) * * 2 / 2 • 

GO TO 10 

20 DO 21 L = 1 f MNMAXO 

S MO=S MO + PH I X ( L ) *#2 
SMT=SMT+PHIT(L)**2 

21 S MF=SMF + PH I (L) **2 

I F < M.GT. MNMAXO) GO TO 11 
SMO=SMO+PHIX(M)**2 
11 BX3(M)=$M0*.5 
BT3(M)=SMT*.5 
BE3 ( M ) = SMF** 5 
BXT3 ( M ) =0« 

GO TO 9 
10 B X3 ( M) = SMO 
BT3 ( M ) = SMT 
BXT3 ( M)=SMR*.5 
BE 3 ( M) = SMF 
9 CONTINUE 
RETURN 
END 


89 



,U3< 


SUBROUTINE EFG(K,MN) 

COMMON 

1/IBL2/NNC 10 ) *MN INIT 
0/BL8/R(200) , GAM <200 ),OMT<200) 

3/BL11/0MXH 200 ) ,PHEE,T0,T2 
4/BL20/DEOMX (200 ) 

4/BL14/LAM2, LSD18,LSD1N 

5/BL15/NU,Ul(10) ,V1< 10) * W1 ( 10) *V2( 10) «U2( 10 ) ,W2(10) 

6 W3 ( 10 ) 

7/BL25/E<4,4) ,F ( 4,4) ,G< 4,4) 

COMMON /BL100/S0RD,TEEO/BL10l/ZO< 4,220) , Z 2 ( 4 , 220) , Z 3 < 4 
1/BL102/DEL0AD/BL103/MASS(200) 

REAL NU»N, LA M2 , LSD1 8 , LSD1 N , MASS , MAS 

N=NN( MN ) 

MAS-MAS S ( K) 

CALL BDB(K,B,DB ,D, DD) 

E(1,1)=B 
E < 1 , 2 ) =0. 

E ( 1 7 3 ) =0* 

E ( 1 7 4 ) =0 • 

E ( 2 , 1 ) =0* 

Dl=< l.-NU) 

RA=R ( K ) 

GA=GAM( K ) 

OX=OMX I ( K) 

OT=OMT( K) 

0 E X=DEOMX ( K ) 

REX= ( 3. *OT-OX) 

G A2=GA**2 
RXE= ( 3. * 0 X— OT) 

OTX-OT *0X 

DNLR=LAM2*D*N*D1/(2.*RA) 

DDNLR=DNLR*DD/D 

E( 2,2)=B*Dl/2.+L AM2*D*Dl*REX**2/8« 

E ( 2 7 3 ) =DNLR*RE X 
E ( 2 7 4) =0* 

E ( 3 7 1 ) =0 • 

E(3,2)=E(2,3) 

R AN= ( N/ RA ) **2 

E(3t3 )=LAM2*D*Dl*(2.*RAN+( l.+NU >*GA2> 

E(3 T 4)=LAM2 
E(4, 1)=0. 

E ( 4 7 2 ) =0 • 

E ( 4 , 3 ) =-D 
E ( 4,4) =0* 

F(ltl)=GA*B+DB 

F< 1,2 ) = (l. + MJ)*B*N/ <2.*RA) + DNLR*REX*RXE/4* 

F ( 1 ,3 ) = B*(OX+NU*OT )+LAM 2 *D*Dl* ( ( 1 *+NU ) *GA 2*0X+RAN*RXE/ 
F ( 1 » 4) =LAM2*0X 
F<2, 1)=-F( 1,2) 

F(2,2)=(Dl/2.) *< GA*B+DB)-(LAM2*D*Dl*REX/8. )*(2.*DEX-GA 
1-3.*0T) )+LAM2*DD*Dl *REX**2/ 8 • 

F( 2,3)=DNLR*(2.*< 1 . +NU ) *GA*OT-DE X+3 .*GA* (OX-OT) ) +DDNLR 
F ( 2 , 4 ) = 0 • 

F ( 3 , 1 ) =-F (1,3) 

F ( 3 , 2 ) =DNLR*( 3 • *GA *OX-GA*OT* (5c+2**NU)-DEX)+DDNLR*REX 
F (3 ,3 )=-LAM2*D*Dl* ( ( 1 .+NU ) *( 2. *GA*OX*OT+G A**3 )+2. *GA*R 
1+LAM2*DD*D1*( ( 1 .+NU ) *GA2+2 .*RAN ) 

F( 3,4>=LAM2*GA*( 2.-NU) 

F ( 4 , 1 ) = D*OX 
F ( 4 , 2 ) = 0 * 

F( 4,3)=-D*NU*GA 
F ( 4 7 4 ) =0* 

G( 1 , 1 )=NU*DB*GA-NU*B*OTX-B*GA2-Dl*B*RAN/2.-LAM2*D*Dl* ( 
10X**2+RXE**2*R AN/8 • ) 

2 - 2.*MAS/DELSD 

G ( 1 , 2 ) =NU*N*DB / RA- ( 3.-NU) / ( 2. *P A) *GA*B*N-DNLR*2.*GA* ( R 
1+(1.+NU)*0TX) 
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G(1»3)=B*(0EX+GA*( OX-OT ) )+DB* (OX+NU*OT)-LAM2*D*Dl*GA*R 
ll.+NU)*OX) 

G ( 1 , 4 ) = L AM2*Dl*GA*OX 

G( 2,1)=-B*GA*N*(3.-NU>/ ( 2.*R A ) - D1*N*DB/ < 2 ,*RA ) +DNLR*2. 
1NU ) *GA*OTX+GA/8 • *( 6.*0TX-7. *OX**2-3. *0T**2 )-DEX/4.*( 5. 
2— DDNLR/4.*REX*RXE 

G( 2, 2 ) =-GA*F ( 2 » 2 ) +D 1/2. *B*OTX-B ! i'R AN-L AM2*0*D1* ( { 1 « + NU) 
l-OTX/8.*REX**2) 

2 - 2.*MAS/DELSD 

G(2»3)=-B*N*( OT+NU*OX )/RA+DNLR*(GA*DEX-2.*GA2#OX-2.*(l 
1*RAN+REX*(GA2+0TX ) ) -DDNLR*REX*GA 
G(2,4)=-NU*LAM2*0T*N/RA 

G( 3» 1) =— B*GA*( OT+NU*OX) +L AM2* D=S=D1 * ( GA#( 1 . +NU ) * ( -G A*DEX 

1- OX*RAN+2.*OTX*OX ) +RAN/2 . *( GA*0X-GA*0T-3. *DEX ) ) 

2- LAM2*DD*Dl*< ( 1 . +NU ) *GA2*OX+RAN/2 . *RXE ) 
G(3,2)=-B*N*(0T+NU*0X) /RA+DNLR*( 2 .*< 1 ,+NU ) * ( OTX*OT-GA2 

l*0T-0T*RAN)+GA*DEX+3. *GA2* ( OT-OX > +OTX*RE X ) -DDNLR* (2.*( 
2*OT+GA*REX) 

G(3,3 )=-B*(OX**2+2. *NU*CTX+OT **2 )+L AM 2*D*D 1#RAN4 ( ( l.+N 
1+2.*GA2 ) +2. * ( G A2 + 0TX ) ) — L AM2*DD=> S D1*R AN* ( 3.+NU ) *GA 
2 - 2«*MAS/DEL SD 

G(3,4)=-LAM2*( D1*0TX+NU*RAN ) 

G(4,1)=D*(DEX+NU*GA*0X) 

G( 4t 2 )= D*NU*N*OT /RA 
G ( 4, 3 ) =D*NU*RAN 
G( 4, 4 ) =— 1. 

RETURN 

END 


SUBROUTINE OUTPUT ( I MODE ) 

REAL NU.MT, MX, MTH , MXT , MTS , KX, KT , K.XT , LAM, LAM2 , MASS 
COMMON /IBL2/N< 1C) ,MNINIT 
2/IBL3/M0,Ml,M2,M3 
1/ I BL4/ KMAX, KL 
2/IBL5/IBCINL,IBCFNL 

3/IBL7/MNMAX0,MAXD( 10) , MAXS( 10) , MAXSYI 10 ) , I S ( 10, 10 ) , JS( 
4,10) ,JD( 10,10) , IJS(IO) 

5/IBL8/LSTEP, ITR 
6/IBL10/I FREQ, NTH MAX 
7/IBL12/KMAX1 ,KMAX2 ,NCONV 
2/IBL13/ITRMAX,LSMAX 

C OM MO N/BL4/P(4, 4,200) , X(4, 200) , ZF 1M ( 4, 4 , 10 ) , ZF2M( 4, 4, 1 
1 ,ZF3M(4,4,10) , ZF4M(4,4, 10) 

COMMON/ BL 5/ TT( 101,MT(10),DT(10) , DMT (10) 

8/BL6/Z( 4,220) , S OE , 0 SE , ALOAD 
9/BL7/D1 , SI 

0/BL8/R ( 200 ) , GAM (200) ,OMT(200) 

3/BL11/0MXI ( 200 ) ,PHEE,T0,T2 
4/BL20/DE0MX( 200 ) 

1/BL10/ PHIX ( 10 ) , PHIT(IO) ,PHI( 10 ) 

3/BL12/TDLI ,TDEL 
4/BL14/L AM2 , LSD18,LSD1N 

5/BL15/NU,Ul( 10 ) ,V1 ( 10) ,W1( 10) , V2< 10) ,U2( 10) ,W2< 10 ) ,U3( 
6W3 ( 10) 

8/BL27/BX3( 10), BT3( 10) ,BXT3( 10) ,BE3( 10) 

COMMON/ BL32/TKN, EL AST, CHAR.SIGO 
1/BL31/DELSQ,EXT1( 10) 

2/BL17/DEL 
3/BL19/TH( 6 ) 

COMMON /BL100/ SORD,TEEO/BL101/Z0(4,220) ,Z2(4,220) ,Z3(4 
1/BL102/DELOAD/BL103/MASS( 200) 

2/BL110/TX(10),TTH( 10) ,TXT( 10) ,MX( 10) , MTH ( 10) ,MXT( 10) ,Q 
3/BLlll/ABZ,ABZ0, AB ZN, ABZ3 , DD2 
DIMENSION PTF(200),PF(200) 

ABZO=SIGO/ELAST 
IF( SORD.NE.O) GO TO 181 
WRITE(6,101 ) LSTEP, ALOAD, ITR 
GO TO 182 

181 TI=LSTEP*DELOAD 
DTI=TI*TEEO 

WRITE (6,151) LSTEP, TI , DTI, ITR 
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CALL POLE ( K ) 
GO TO 999 


182 LAM=TKN/CHAR 
ENL=1 % 

ABZ=SIGO*TKN 
ABZ3 =ABZ*TKN*TKN/CHAR 
ABZN=CHAR*SIGO/ELAST 
IF ( ITRMAX. E Q* 1 ) ENL=0. 

DD2=1 •- NU*^2 
D2 I =1 • / DD2 
DP 1= 1. / S 1 
DNI=1./D1 
TDLSQI=.5/DELSQ 
I F ( NTHM AX. EQ. 0 ) GO TO 991 
DO 21NTH=1 » NTHMAX 
DO 1 MN = 1 » MNMAXO 
Il=l+< MN-1)*KMAX2 
12=11+1 

U1(MN)-Z(ly ID 
U2 ( MN ) = Z ( 1 1 12) 

VI ( MN ) =Z ( 2 » 1 1 ) 

V2(MIM) = Z(2» 12) 

W1(MN)=Z(3, II ) 

1 W2 ( MN ) = Z ( 3 t 12) 

THET=TH( NTH) 

WRITE(6,116) THET 
DO 121 K=1,KMAX 
K 1 = K+1 

CALL BDB(K,BS,DB,DS,DD> 

IF(K.EQ. 1. AND. IBCINL.LT. 0) CALL POLE(K) 

IF<K. EQ.l. AND. IBCINL.LT. 0) GO TO 999 
IF( K.EQ.KMAX.AND. IBCFNL.LT.O) 

IF( K.EQ. KM AX. AND. I BCFNL.LT. 0) 

CALL PH I BET ( K ) 

DEX=DEOMX(K) 

RRA=1./R(K) 

OX=OMX I ( K ) 

OT=OMT(K) 

GA=GAM ( K ) 

DCXT=OX-OT 
GDO=GA*DOXT 
DD2D=DD 2*DS 
DO 3 MN = 1 y MNMAXO 
EN=N( MN ) 

ENR=EN*RRA 
CALL TLOAD(K) 

TTS=TT ( MN) * ALOAD 

EX=(U3(MN)-U1(MN) )*TDLI +0X*W2<MN) + ENL*OSE* ( BX3 ( MN ) + 
ET=ENR# V2 ( MN ) + GA*U2(MN) + 0T*W2(MN) + ENL^OSE’M BT3 < M 
EXT=.5*( (V3(MN)-V1( MN) )*TDLI - FNR*U2(MN) - GA*V2<MN) 

1 + ENL*S0E*BXT3(MN) ) 

KT=ENR#PHIT (MN) + GA*PHIX(MN) 

KXT=.5*(ENR*(-PHIX(MN) -GA*W2(MN) + ( W3 < MN ) -W 1 ( MN ) ) *T D 

1 + GD0*V2(MN) + OT* ( V3 ( MN ) -V 1 ( MN ) ) *TDL I 

2 -GA*PHIT(MN)-DOXT*PHI (MN) ) 

TXt MN ) = BS* ( E X+NU*E T )-TTS 

TTH( MN ) = BS*( ET+NU*EX J-TTS 
T XT ( MN )=BS*D1*EXT 
MK1=K1+ ( MN- 1 ) *KMAX2 
MX ( MN ) = Z ( 4 * MK1 ) 

MTH( MN) =NU*MX( MN ) +DD2D*KT-D1*MT ( MN ) 5 
MXT < MN)=DS*D1*KXT 
MK1 1=MK 1+1 
MKK 1=MK 1—1 

QS ( MN ) =S IG0*TKN*LAM2*( GA*MX ( MN ) +( Z ( 4« MK11 )-Z(4,MKKl ) )* 
1 +ENR*MXT ( MN ) -GA*MTH ( MN ) ) 

MX(MN)=MX(MN)*ABZ3 
MTH(MN) = MTH ( MN ) * ABZ3 
MXT ( MN)=MXT(MN)*ABZ3 
TX(MN)=TX(MN)*ABZ 
TTH( MN ) =TTH( MN ) *ABZ 
TXT(MN)=TXT(MN)*ABZ 
PH I X ( MN )=PHIX ( MN ) *ABZO 
PH I T ( MN )=PHIT ( MN ) * ABZO 


^ ALOAD 
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PHI (MN)=PHI <MN)*ABZO 
U1(MN)=U2(MN) 

U2 ( MN ) =U3 ( MN) 

Vl(MN) =V2 ( MN ) 

V2 ( MN ) =V3 ( MN) 

W1 ( MN ) =W2 ( MN) 

3 W2 ( MN ) =W3 ( MN ) 

FK=K— X 

FI FREQ=I FREQ 

KTST=(K>1)/IFREQ 

FKTST=KT$T 

FKTEST=FK/FI FREQ-FKTST 
IF(K.EQ.l.OR.K.EQ.KMAX) GO TO 999 
I F ( FKTEST. NE.O • ) GO TO 2 

999 


PTF ( K ) =0 • 

PF ( K ) =0 • 

AMX=0. 

AMTH=0 • 

AMXTH=0. 

ANX=0. 

ANTH=0 • 

ANXTH=0 • 

AQS = 0.0 
DO 72 MN = 1 y MNM AXO 
EN=N( MN ) 

FC=EN*T HET 
SN=SI N ( FC ) 

C S=COS ( FC ) 

X( 1 ,K)=X(l t K)+U1 (MN )*CS*ABZN 
X(2fK)=X(2»K)+Vl (MN)*SN*ABZN 
X( 3,K)=X(3,K)+W1(MN)*CS*ABZN 
X(4,K)=X(4,K)+PHIX(MN)*CS 
PTF ( K ) =PTF ( K ) +PH IT ( MN) *SN 
AMX=AMX+MX(MN)*CS 
AMTH=AMTH+MTH( MN)*CS 
AMXTH=AMXTH+MXT ( MN) *SN 
ANX=ANX+TX(MN)*CS 
ANTH=ANTH+TTH( MN)*CS 

AQS = AQS + QS ( MN ) *C S 
ANXTH=A NXTH+TXT (MN)*SN 
72 PF(K)=PF(K)+PHI <MN)*SN 
429 CONTINUE 

IF(K.EO.l) W RI T E ( 6 » 117) 

117 FORMAT ( / » 8H ST ATI0N15H N S 16H N THETA 

IN STHETA 16H Q S 16H MS 16H 

216H M STHETA // ) 

WRITE ( 6, 118) KtANX, ANTH,ANXTH , AQS , AMX , A MTH , AMXTH 

118 FORMAT ( lXt I3f3X t 7E16.4) 

101 FORMAT ( 1 HI » 34H THE LOAD STEP NUMBER IS I2,29H 

IE LOAD FACTOR IS E11.4,36H THE SOLUTION CONVE 

2 1H ITERATIONS////) 

151 FORMAT ( 1H1 ,26H THE TIME STEP NUMBER IS I4,18H TH 

1.2»4H OR E9.3,40H SECONDS THE SOLUTION CONVERGED 

3ERAT IONS/// / ) 

102 FORMAT ( 1H0 1 35H MODAL OUTPUT FOR STATION 13//) 

103 FORMAT (//7H N 16H U 16H V 

1W 16H PHI S 16H PHI THETA10H 

104 FORMAT ( 2X»I2?3X»6E16.4) 

105 FORMAT ( 7H N 16H N S 16H N THETA 

1 STHETA 16H Q S 16H MS 16H 

2H M STHETA /) 

106 FORMAT ( 1 HO * 8HS T AT I ON I3,14H- THETA OUTPUT/ ) 

107 FORMAT ( 14H THETA 14H U 

1 V 14H W 14H PHI S 14H PH 

2 PHI ) 

108 FORMAT ( 8E14. 4) 

109 FORMAT (lHOt 14H THETA 14H M 

1 M THETA 14H M STHETA 14H N S 14H 


X( 1 « K ) =0. 
X ( 2 ? K ) =0 • 
X(3,K)=0. 
X(4tK)=0. 
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2H N STHETA 

114 FORMAT ( 2X, 12, 3X 

14H Q 

, 7E 1 6. 4 ) 

S 

) 


115 FORMAT ( 7H N 

16H 

N S 

1 6 H 

N THETA 

1STHETA 16H 

M S 

16H 

M THETA 

15H 


C 116 FORMAT ( 1H1 , 84H THE SUMMED FORCES, MOMENTS, DI 

116 FORMAT ( 1HO ,84H THE SUMMED FORCES, MOMENTS, DI 

1 AND ROTATIONS FOLLOW FOR THETA =E15.6///) 

2 CONTINUE 
121 CONTINUE 

DO 660 K=1 , KMAX 
FK=K-1 

F IFREQ = I FREQ 

KTST=(K-1)/IFREQ 

FKTST=KTST 

FKTEST=FK/F IFREQ-FKTST 
IF(K.EQ.l.OR.K.EQ.KMAX) GO TO 661 
IFCFKTEST.NE.O. ) GO TO 658 
661 IF(K.EQ.l) WRI TE (6,217) 

WRITE (6,218) K,X(1 ,K) ,X(2,K) ,X(3,K) ,X(4,K),PTF(K) ,PF(K 

217 FORMAT ( //8H STATION, 15H U 16H V 

1 W 16H PHI S 16H PHI THETA16H 

2 /) 

218 FORMAT ( 1X,I3,3X,6E16.4) 

658 DO 659 1=1,4 

659 X( I , K) =0# 

660 CONTINUE 
21 CONTINUE 

991 IF< IMODE.LE.O) RETURN 
DO 534 MN=1 ,MNMAXO 
WRITE(6,749) N(MN) 

749 FORMAT( 1H1,40X,27H MODAL OUTPUT FOR MODE N = 13, 8H FOL 
DO 521 MM=1 , MNMAXO 
I l = l + ( MM— 1 ) *KM AX2 
12=11+1 

U1(MM) = Z( 1, ID 
U2 ( MM ) =Z ( 1 , 12 ) 

VI ( MM ) = Z ( 2 , II ) 

V2(MM) = Z(2, 12) 

W1 ( MM ) = Z ( 3 , 1 1 ) 

W2(MM)=Z(3, 12) 

521 CONTINUE 

DO 445 K=1 , KMAX 
K 1=K+1 

CALL BDB(K,BS,DB,DS,DD) 

IF(K.EQ. l.AND. IBCINL.LT.O) CALL POLE(K) 

IF(K.EQ. KMAX. AND. IBCFNL.LT.O) CALL POLE(K) 

TXZ=TX ( MN ) 

TTHZ=TTH(MN) 

TXTZ=TXT ( MN ) 

AMXZ=MX ( MN ) 

AMTHZ=MTH( MN ) 

AMXTZ=MXT ( MN ) 

QSZ=QS ( MN) 

X ( 1 , K ) =PHI X ( MN ) 

X ( 2 , K ) = PHIT ( MN ) 

X ( 3 , K) = PHI ( MN) 

IF(K.EQ.l .AND. IBCINL.LT.O) GO TO 583 
IF(K.EQ. KMAX. AND. IBCFNL.LT.O) GO TO 583 
CALL PH I BET (K) 

DEX=DEOMX( K ) 

RRA=1 • / R ( K ) 

OX=OMXI ( K) 

OT=OMT(K) 

GA=GAM ( K ) 

DOXT = 0X— OT 
GD0=GA*D0XT 
DD2D=DD2*DS 
EN=N( MN ) 

ENR=EN*RRA 
CALL TLOAD(K) 

TTS=TT(MN)*ALOAD 

EX=(U3(MN)-U1(MN))*TDLI +0X*W2(MN) + ENL *0SE* ( B X3 ( MN ) + 
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533 


583 


445 


593 


446 

534 


ET=ENR*V2(MN) + GA*U2(MN) + 0T*W2(MN) + ENL*OSE* ( BT3 ( M 
EXT=.5*( <V3(MN)-V1< MN) )*TDLI - ENR*U2 (MN) - GA*V2(MN) 

1 + ENL*S0E*BXT3 ( MN ) ) 

KT =ENR* PHI T ( MN) + GA*PHIX(MN) 

KXT=.5*< ENR*(-PHIX(MN) -GA*W2(MN) + ( W3 ( MN ) -W1 ( MN ) ) *TD 

1 + GD0*V2 ( MN ) + 0T*(V3(MN)-V1 (MN))*TDII 

2 — GA*PHIT(MNJ — DOXT*PH I ( MN ) ) 

TXZ = <BS*(EX+NU*ET)-TTS)*ABZ 
TTHZ = ( BS*(ET+NU*EX)-TTS)*ABZ 
TXTZ = BS*D1*EXT*ABZ 
MK1=K1+ ( MN-1 ) *KMAX2 

AMXZ = Z ( 4» MK1 ) 

AMTHZ = NU*AMXZ+DD2D*KT-D1*MT(MN)*AL0AD 

AMXTZ = DS*D1* KXT 

MK11=MK1+1 

LI 1/ 1/ *J ^ |J 1/ 1 _ 1 

QSZ = SIG0*TKN*LAM2*(GA*AMXZ +( Z ( 4 t MK1 1 ) -Z ( 4, MKKl ) >* 
1 +ENR*AMXTZ -GA*AMTHZ) 

AMXZ=AMXZ*ABZ3 
AMTHZ= AMTHZ*AB Z3 
AMXTZ=AMXTZ*ABZ3 
X(1,K) = PHIX (MN)*ABZO 
X ( 2 » K ) = PH IT ( MN ) *AB ZO 
X ( 3 » K ) *PHI ( MN ) * ABZO 
DO 533 MM=1 f MNMAXO 
U 1 < MM ) =U2 ( MM) 

U2 ( MM ) =U3 ( MM ) 

Vl(MM) =V2( MM) 

V2 ( MM ) =V 3 ( MM ) 

Wl(MM) =W2(MM) 

W2(MM)=W3(MM) 

FK=K-1 

F I FREQ= I FREQ 
KTST=(K— 1)/IFREQ 
FKTST=KTST 

FKTEST=FK/FI FREQ-FKTST 
IF(K«EQ*l#OR*K*EQ#KMAX) GO TO 583 
I F C FKTEST.NE.O. ) GO TO 445 
CONTINUE 

IF(K.EQ.l) WRITE(6t 117) 

WRITE (6, 118 ) K, TXZ T TTHZ, TXTZ ,QS Z , AMXZ t AMT HZ , AMXTZ 

CONTINUE 

WRITE<6,217 ) 

DO 446 K=1,KMAX 
FK=K— 1 

F I FREQ= I FREQ 
KTST= ( K-l ) /IFREQ 
FKTST = KT ST 

FKTEST=FK/FIFRE Q-FKTST 

IF(K.EQ. l.OR.K. EQ.KMAX) GO TO 593 

IF(FKTEST*NE#0* ) GO TO 446 

KZ=K+1+ ( MN— 1 ) #KMAX2 

UP = Z( l t KZ)*ABZN 

VP=Z(2, KZ)*ABZN 

WP=Z(3,KZ)*ABZN 

WRITE (6*218) K*UP, VP*WP*X<1 *K) ,X(2,K) ,X(3,K) 

CONTINUE 

CONTINUE 

RETURN 

END 
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TABLE I. IMPORTANT FORTRAN VARIABLES 


FORTRAN Variable 

A(4,4) 

bee(4,4) 

BEl (10) 

BE2(lO) 

BE3(10) 

BTl(lO) 

BT2 (10) 

BT3(10) 

BXl(lO) 

BX2 (10) 

BX3(10) 

BXT1(10) 

BXT2(10) 

BXT3(10) 

c(4,4) 

capll(4,4) 

cafli(4,4) 

DEE(4,4,200) 

DELSD 

DST(4, 4,200) 

D1 

E(4,4) 

ell(4) 

eli(4) 

ETl(lO) 

ET2 (10) 

ET3(10) 

ETTl(lO) 

ETT2 (10) 

ETT3(10) 


Definition 

A matrix 

B matrix 

p at i - 1, 
i, and i + 1 

Pg at i - 1, 

i, and i + 1 

P at i - 1. 
s 5 

i , _and i + 1 

P sQ at i - 1, 
i 5 and i + 1 

C_ matrix 

Ajj, A k 

Al’ Al 

["b. - C.P. ,1 

L i i l-l.l 
DEL0AD*DEL0AD 

[ s i - 

1 - V 

E matrix 

Al 

TIq at i - 1, 
1 , and i + 1 

n ee at i - 1 
i, and i + 1 

4 0s at i - 1 
i, and i + 1 


Em (io) 

ETX2(10) 

ETX3(10) 
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FORTRAN Variable 


Definition 


EX1(10) 

T| s at i - 1, 

EX2(10) 

EX3(10) 

1, and i + 1 

EXT1(10) 

Tl se at i - 1 

EXT2(10) 

EXT3(10) 

i, and i + 1 

EXXl(lO) 

V at 1 - 1, 

EXX2(lO) 

EXX3(10) 

i 5 and i + 1 

F(4,4) 

F matrix 

FFS(4, 10) 

f^ matrix 

fls(4) 

f matrix 

G(4,4) 

G matrix 

gee(4) 

g matrix 

GEES (4,10) 

matrix 

H(4,4) 

H matrix 

ID ( 10 , 10 ) 

See description of 
subroutine MODES 

IJS(IO) 

See description of 
subroutine MODES 

IS (10,10) 

See description of 
subroutine MODES 

ITR 

Iteration number 

jay(4,4) 

J matrix 

jd(io,io) 

See description of 
subroutine MODES 
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FORTRAN Variable 


Definition 


KL 

KMAX-1 

KLL 

KL-1 

KMAX1 

KMAX+1 

KMAX2 

KMAX+2 

js(io,io) 

See description of 
subrountine MODES 

MAXD(lO) 

See description of 
subroutine MODES 

MAXS(lO) 

See description of 
subroutine MODES 

MAXSY(lO) 

See description of 
subroutine MODES 

MO, Ml, V2., and M3 

r\ 

H ro 

li 11 

ss 

ri fN 

O OJ 
II II 

11 

N(10) or NN(10) 

n 

0MEGL(4,4) 

°K 

0megi(4,4) 


0SE 

.5° /E 
o' 0 

P(4,4,200) 

P matrix 

PHI (10) 

c p 

PHIT(IO) 

CD 

-Q 

CD 

PHIX(lO) 

$ j Cp 

s 

S0E 

cr /E 
o' 0 

SI 

1 + V 
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FORTRAN Variable 

Definition 

TDEL 

2A 

TDLI 

1/(2A) 

th(6) 

e 

Ul(lO) 

U, u at i - ] 

U2(l0) 

i 5 and i + 1 

U3(10) 


Vl(lO) 


V2(l0 

V, v at i ■ ] 

V3(10) 

± 9 and i + 1 

W1(10) 

w at i - ] 

W2(10) 

i* and i + 1 

W3(10) 


X(4,200) 

x matrix 

Z(4,220) 

z matrix 

ZD0T(4,22O) 

dz/dt 

Z0(4,22O) 

z at t - fit 

Z2(4,220) 

z at t - 2 fit 

Z3(4,220) 

z at t - 36t 


I 




Figure 10. Flow of Program Logic in MAIN 
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